#This is the replication file for Thomann, E. van Engen, N. and L. Tummers. Forthcoming. The necessity of discretion: a behavioral evaluation of bottom-up implementation theory.  Journal of Public Administration Research and Theory.

# This code replicates the analysis of dataset 1, direct calibration method.

#R code based  on Thomann, Eva and Stefan Wittwer (2016). Performing fuzzy- and crisp set QCA with R: A user-oriented beginner's guide. Version 19.9.2016. URL: http://www.evathomann.com/links/qca-r-manual [date of access: 25.11.2016].


#clear workspace
rm(list = ls())


library(lattice); library(arm); library(xtable); library(foreign); library(psych); library(directlabels); library(betareg); library(VIM); library(base); library(plyr); library(dplyr); library(QCA); library(QCAGUI); library(SetMethods)

setwd("C:/Users/Eva/Dropbox/Arbeit/Publikationen/Papers/Two-factor theory/Material/Data/Dataset 1")


##############PREPARING THE DATASET###########

tfr <- read.csv("tfr_d1.csv", row.names=1, header = TRUE, sep = ";",  dec = ",")
describe(tfr)

#check missings
missings <- as.numeric(complete.cases(tfr))
tabmiss <- table(missings)
prop.table(tabmiss)

#identify  cases with missings on all indicators for a causal factor
##for W
colnames(tfr)

tfr$w1miss <- as.numeric(is.na(tfr$W1))
tabw1miss <- table(tfr$w1miss)
prop.table(tabw1miss)


tfr$w2miss <- as.numeric(is.na(tfr$W2))
tabw2miss <- table(tfr$w2miss)
prop.table(tabw2miss)

tfr$w3miss <- as.numeric(is.na(tfr$W3))
tabw3miss <- table(tfr$w3miss)
prop.table(tabw3miss)

tfr$w4miss <- as.numeric(is.na(tfr$W4))
tabw4miss <- table(tfr$w4miss)
prop.table(tabw4miss)

tfr$w5miss <- as.numeric(is.na(tfr$W5))
tabw5miss <- table(tfr$w5miss)
prop.table(tabw5miss)

colnames(tfr)
tfr$wmiss <- rowSums(tfr[,50:54])
tfr$wmiss


tfr$wmissout <- as.numeric(tfr$wmiss ==5)
tfr$wmissout
tabwmissout <- table(tfr$wmissout)
prop.table(tabwmissout)

colnames(tfr)

##for SP

tfr$sp1miss <- as.numeric(is.na(tfr$SP1))
tabsp1miss <- table(tfr$sp1miss)
prop.table(tabsp1miss)


tfr$sp2miss <- as.numeric(is.na(tfr$SP2))
tabsp2miss <- table(tfr$sp2miss)
prop.table(tabsp2miss)

tfr$sp3miss <- as.numeric(is.na(tfr$SP3))
tabsp3miss <- table(tfr$sp3miss)
prop.table(tabsp3miss)

tfr$sp4miss <- as.numeric(is.na(tfr$SP4))
tabsp4miss <- table(tfr$sp4miss)
prop.table(tabsp4miss)

tfr$sp5miss <- as.numeric(is.na(tfr$SP5))
tabsp5miss <- table(tfr$sp5miss)
prop.table(tabsp5miss)

tfr$sp6miss <- as.numeric(is.na(tfr$SP6))
tabsp6miss <- table(tfr$sp6miss)
prop.table(tabsp6miss)

colnames(tfr)

tfr$spmiss <- rowSums(tfr[,57:62])
tfr$spmiss

tfr$spmissout <- as.numeric(tfr$spmiss==6)
tfr$spmissout
tabspmissout <- table(tfr$spmissout)
prop.table(tabspmissout)

##for TP
tfr$tp1miss <- as.numeric(is.na(tfr$TP1))
tabtp1miss <- table(tfr$tp1miss)
prop.table(tabtp1miss)


tfr$tp2miss <- as.numeric(is.na(tfr$TP2))
tabtp2miss <- table(tfr$tp2miss)
prop.table(tabtp2miss)

tfr$tp3miss <- as.numeric(is.na(tfr$TP3))
tabtp3miss <- table(tfr$tp3miss)
prop.table(tabtp3miss)

tfr$tp4miss <- as.numeric(is.na(tfr$TP4))
tabtp4miss <- table(tfr$tp4miss)
prop.table(tabtp4miss)

tfr$tp5miss <- as.numeric(is.na(tfr$TP5))
tabtp5miss <- table(tfr$tp5miss)
prop.table(tabtp5miss)

tfr$tp6miss <- as.numeric(is.na(tfr$TP6))
tabtp6miss <- table(tfr$tp6miss)
prop.table(tabtp6miss)

colnames(tfr)

tfr$tpmiss <- rowSums(tfr[,65:70])
tfr$tpmiss

tfr$tpmissout <- as.numeric(tfr$tpmiss==6)
tfr$tpmissout
tabtpmissout <- table(tfr$tpmissout)
prop.table(tabtpmissout)

##for OP
tfr$op1miss <- as.numeric(is.na(tfr$OP1))
tabop1miss <- table(tfr$op1miss)
prop.table(tabop1miss)

tfr$op2miss <- as.numeric(is.na(tfr$OP2))
tabop2miss <- table(tfr$op2miss)
prop.table(tabop2miss)

tfr$op3miss <- as.numeric(is.na(tfr$OP3))
tabop3miss <- table(tfr$op3miss)
prop.table(tabop3miss)

tfr$op4miss <- as.numeric(is.na(tfr$OP4))
tabop4miss <- table(tfr$op4miss)
prop.table(tabop4miss)

tfr$op5miss <- as.numeric(is.na(tfr$OP5))
tabop5miss <- table(tfr$op5miss)
prop.table(tabop5miss)

tfr$op6miss <- as.numeric(is.na(tfr$OP6))
tabop6miss <- table(tfr$op6miss)
prop.table(tabop6miss)

colnames(tfr)

tfr$opmiss <- rowSums(tfr[,73:78])
tfr$opmiss

tfr$opmissout <- as.numeric(tfr$opmiss==6)
tfr$opmissout
tabopmissout <- table(tfr$opmissout)
prop.table(tabopmissout)

##for SM

tfr$sm1miss <- as.numeric(is.na(tfr$SM1))
tabsm1miss <- table(tfr$sm1miss)
prop.table(tabsm1miss)

tfr$sm2miss <- as.numeric(is.na(tfr$SM2))
tabsm2miss <- table(tfr$sm2miss)
prop.table(tabsm2miss)

tfr$sm3miss <- as.numeric(is.na(tfr$SM3))
tabsm3miss <- table(tfr$sm3miss)
prop.table(tabsm3miss)

tfr$sm4miss <- as.numeric(is.na(tfr$SM4))
tabsm4miss <- table(tfr$sm4miss)
prop.table(tabsm4miss)

tfr$sm5miss <- as.numeric(is.na(tfr$SM5))
tabsm5miss <- table(tfr$sm5miss)
prop.table(tabsm5miss)

tfr$sm6miss <- as.numeric(is.na(tfr$SM6))
tabsm6miss <- table(tfr$sm6miss)
prop.table(tabsm6miss)

tfr$sm7miss <- as.numeric(is.na(tfr$SM7))
tabsm7miss <- table(tfr$sm7miss)
prop.table(tabsm7miss)

tfr$sm8miss <- as.numeric(is.na(tfr$SM8))
tabsm8miss <- table(tfr$sm8miss)
prop.table(tabsm8miss)

tfr$sm9miss <- as.numeric(is.na(tfr$SM9))
tabsm9miss <- table(tfr$sm9miss)
prop.table(tabsm9miss)

tfr$sm10miss <- as.numeric(is.na(tfr$SM10))
tabsm10miss <- table(tfr$sm10miss)
prop.table(tabsm10miss)

tfr$sm11miss <- as.numeric(is.na(tfr$SM11))
tabsm11miss <- table(tfr$sm11miss)
prop.table(tabsm11miss)

tfr$sm12miss <- as.numeric(is.na(tfr$SM12))
tabsm12miss <- table(tfr$sm12miss)
prop.table(tabsm12miss)

colnames(tfr)

tfr$smmiss <- rowSums(tfr[,81:92])
tfr$smmiss

tfr$smmissout <- as.numeric(tfr$smmiss==12)
tfr$smmissout
tabsmmissout <- table(tfr$smmissout)
prop.table(tabsmmissout)

##for CM

tfr$cm1miss <- as.numeric(is.na(tfr$CM1))
tabcm1miss <- table(tfr$cm1miss)
prop.table(tabcm1miss)


tfr$cm2miss <- as.numeric(is.na(tfr$CM2))
tabcm2miss <- table(tfr$cm2miss)
prop.table(tabcm2miss)

tfr$cm3miss <- as.numeric(is.na(tfr$CM3))
tabcm3miss <- table(tfr$cm3miss)
prop.table(tabcm3miss)

tfr$cm4miss <- as.numeric(is.na(tfr$CM4))
tabcm4miss <- table(tfr$cm4miss)
prop.table(tabcm4miss)

colnames(tfr)

tfr$cmmiss <- rowSums(tfr[,95:98])
tfr$cmmiss

tfr$cmmissout <- as.numeric(tfr$cmmiss ==4)
tfr$cmmissout
tabcmmissout <- table(tfr$cmmissout)
prop.table(tabcmmissout)

#delete cases with missings on all indicators for 1 causal factor

tfr = tfr[tfr$wmissout == 0,]
tfr = tfr[tfr$spmissout == 0,]
tfr = tfr[tfr$tpmissout == 0,]
tfr = tfr[tfr$opmissout == 0,]
tfr = tfr[tfr$smmissout == 0,]
tfr = tfr[tfr$cmmissout == 0,]
tfr

describe(tfr)

#check missings
missings <- as.numeric(complete.cases(tfr))
tabmiss <- table(missings)
prop.table(tabmiss)

######Calibration############
#calibrate indicator sets, recode neutral answers on the Likert scale into full nonmembership; calibrate indicator sets
#the values on the likert scales for powerlessness and meaninglessness need to be reverted because they should mean 
#powerfulness and meaningfulness
##direct method
###for W


tfr$wa <- tfr$W1
tfr$wa[tfr$W1 == 4] <- 3
tfr$wa[tfr$W1 == 5] <- 4
tfr$wa


tfr$wb <- tfr$W2
tfr$wb[tfr$W2 == 4] <- 3
tfr$wb[tfr$W2 == 5] <- 4

tfr$wc <- tfr$W3
tfr$wc[tfr$W3 == 4] <- 3
tfr$wc[tfr$W3 == 5] <- 4


tfr$wd <- tfr$W4
tfr$wd[tfr$W4 == 4] <- 3
tfr$wd[tfr$W4 == 5] <- 4


tfr$we <- tfr$W5
tfr$we[tfr$W5 == 4] <- 3
tfr$we[tfr$W5 == 5] <- 4


tfr$w1 <- calibrate(tfr$wa, type = "fuzzy", thresholds = c(1, 2.5, 4), logistic = TRUE)
tfr$w1[tfr$W1 == 3] <- 0.05
tfr$w1
tfr$w2 <- calibrate(tfr$wb, type = "fuzzy", thresholds = c(1, 2.5, 4), logistic = TRUE)
tfr$w2[tfr$W2 == 3] <- 0.05
tfr$w2
tfr$w3 <- calibrate(tfr$wc, type = "fuzzy", thresholds = c(1, 2.5, 4), logistic = TRUE)
tfr$w3[tfr$W3 == 3] <- 0.05
tfr$w3
tfr$w4 <- calibrate(tfr$wd, type = "fuzzy", thresholds = c(1, 2.5, 4), logistic = TRUE)
tfr$w4[tfr$W4 == 3] <- 0.05
tfr$w4
tfr$w5 <- calibrate(tfr$we, type = "fuzzy", thresholds = c(1, 2.5, 4), logistic = TRUE)
tfr$w5[tfr$W5 == 3] <- 0.05
tfr$w5

## for SP
tfr$spa <- tfr$SP1
tfr$spa[tfr$SP1 == 4] <- 3
tfr$spa[tfr$SP1 == 5] <- 4
tfr$spa


tfr$spb <- tfr$SP2
tfr$spb[tfr$SP2 == 4] <- 3
tfr$spb[tfr$SP2 == 5] <- 4

tfr$spc <- tfr$SP3
tfr$spc[tfr$SP3 == 4] <- 3
tfr$spc[tfr$SP3 == 5] <- 4


tfr$spd <- tfr$SP4
tfr$spd[tfr$SP4 == 4] <- 3
tfr$spd[tfr$SP4 == 5] <- 4


tfr$spe <- tfr$SP5
tfr$spe[tfr$SP5 == 4] <- 3
tfr$spe[tfr$SP5 == 5] <- 4

tfr$spf <- tfr$SP6
tfr$spf[tfr$SP6 == 4] <- 3
tfr$spf[tfr$SP6 == 5] <- 4


tfr$sp1 <- calibrate(tfr$spa, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sp1[tfr$SP1 == 3] <- 0.05
tfr$sp1
tfr$sp2 <- calibrate(tfr$spb, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sp2[tfr$SP2 == 3] <- 0.05
tfr$sp2
tfr$sp3 <- calibrate(tfr$spc, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sp3[tfr$SP3 == 3] <- 0.05
tfr$sp3
tfr$sp4 <- calibrate(tfr$spd, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sp4[tfr$SP4 == 3] <- 0.05
tfr$sp4
tfr$sp5 <- calibrate(tfr$spe, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sp5[tfr$SP5 == 3] <- 0.05
tfr$sp5
tfr$sp6 <- calibrate(tfr$spf, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sp6[tfr$SP6 == 3] <- 0.05
tfr$sp6

## for TP
tfr$tpa <- tfr$TP1
tfr$tpa[tfr$TP1 == 4] <- 3
tfr$tpa[tfr$TP1 == 5] <- 4
tfr$tpa


tfr$tpb <- tfr$TP2
tfr$tpb[tfr$TP2 == 4] <- 3
tfr$tpb[tfr$TP2 == 5] <- 4

tfr$tpc <- tfr$TP3
tfr$tpc[tfr$TP3 == 4] <- 3
tfr$tpc[tfr$TP3 == 5] <- 4


tfr$tpd <- tfr$TP4
tfr$tpd[tfr$TP4 == 4] <- 3
tfr$tpd[tfr$TP4 == 5] <- 4


tfr$tpe <- tfr$TP5
tfr$tpe[tfr$TP5 == 4] <- 3
tfr$tpe[tfr$TP5 == 5] <- 4

tfr$tpf <- tfr$TP6
tfr$tpf[tfr$TP6 == 4] <- 3
tfr$tpf[tfr$TP6 == 5] <- 4


tfr$tp1 <- calibrate(tfr$tpa, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$tp1[tfr$TP1 == 3] <- 0.05
tfr$tp1
tfr$tp2 <- calibrate(tfr$tpb, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$tp2[tfr$TP2 == 3] <- 0.05
tfr$tp2
tfr$tp3 <- calibrate(tfr$tpc, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$tp3[tfr$TP3 == 3] <- 0.05
tfr$tp3
tfr$tp4 <- calibrate(tfr$tpd, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$tp4[tfr$TP4 == 3] <- 0.05
tfr$tp4
tfr$tp5 <- calibrate(tfr$tpe, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$tp5[tfr$TP5 == 3] <- 0.05
tfr$tp5
tfr$tp6 <- calibrate(tfr$tpf, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$tp6[tfr$TP6 == 3] <- 0.05
tfr$tp6

##for OP

tfr$opa <- tfr$OP1
tfr$opa[tfr$OP1 == 4] <- 3
tfr$opa[tfr$OP1 == 5] <- 4
tfr$opa


tfr$opb <- tfr$OP2
tfr$opb[tfr$OP2 == 4] <- 3
tfr$opb[tfr$OP2 == 5] <- 4

tfr$opc <- tfr$OP3
tfr$opc[tfr$OP3 == 4] <- 3
tfr$opc[tfr$OP3 == 5] <- 4


tfr$opd <- tfr$OP4
tfr$opd[tfr$OP4 == 4] <- 3
tfr$opd[tfr$OP4 == 5] <- 4


tfr$ope <- tfr$OP5
tfr$ope[tfr$OP5 == 4] <- 3
tfr$ope[tfr$OP5 == 5] <- 4

tfr$opf <- tfr$OP6
tfr$opf[tfr$OP6 == 4] <- 3
tfr$opf[tfr$OP6 == 5] <- 4


tfr$op1 <- calibrate(tfr$opa, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$op1[tfr$OP1 == 3] <- 0.05
tfr$op1
tfr$op2 <- calibrate(tfr$opb, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$op2[tfr$OP2 == 3] <- 0.05
tfr$op2
tfr$op3 <- calibrate(tfr$opc, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$op3[tfr$OP3 == 3] <- 0.05
tfr$op3
tfr$op4 <- calibrate(tfr$opd, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$op4[tfr$OP4 == 3] <- 0.05
tfr$op4
tfr$op5 <- calibrate(tfr$ope, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$op5[tfr$OP5 == 3] <- 0.05
tfr$op5
tfr$op6 <- calibrate(tfr$opf, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$op6[tfr$OP6 == 3] <- 0.05
tfr$op6

##for SM

tfr$sma <- tfr$SM1
tfr$sma[tfr$SM1 == 4] <- 3
tfr$sma[tfr$SM1 == 5] <- 4
tfr$sma


tfr$smb <- tfr$SM2
tfr$smb[tfr$SM2 == 4] <- 3
tfr$smb[tfr$SM2 == 5] <- 4

tfr$smc <- tfr$SM3
tfr$smc[tfr$SM3 == 4] <- 3
tfr$smc[tfr$SM3 == 5] <- 4


tfr$smd <- tfr$SM4
tfr$smd[tfr$SM4 == 4] <- 3
tfr$smd[tfr$SM4 == 5] <- 4


tfr$sme <- tfr$SM5
tfr$sme[tfr$SM5 == 4] <- 3
tfr$sme[tfr$SM5 == 5] <- 4

tfr$smf <- tfr$SM6
tfr$smf[tfr$SM6 == 4] <- 3
tfr$smf[tfr$SM6 == 5] <- 4


tfr$smg <- tfr$SM7
tfr$smg[tfr$SM7 == 4] <- 3
tfr$smg[tfr$SM7 == 5] <- 4
tfr$smg


tfr$smh <- tfr$SM8
tfr$smh[tfr$SM8 == 4] <- 3
tfr$smh[tfr$SM8 == 5] <- 4

tfr$smi <- tfr$SM9
tfr$smi[tfr$SM9 == 4] <- 3
tfr$smi[tfr$SM9 == 5] <- 4


tfr$smj <- tfr$SM10
tfr$smj[tfr$SM10 == 4] <- 3
tfr$smj[tfr$SM10 == 5] <- 4


tfr$smk <- tfr$SM11
tfr$smk[tfr$SM11 == 4] <- 3
tfr$smk[tfr$SM11 == 5] <- 4

tfr$sml <- tfr$SM12
tfr$sml[tfr$SM12 == 4] <- 3
tfr$sml[tfr$SM12 == 5] <- 4

tfr$sm1 <- calibrate(tfr$sma, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sm1[tfr$SM1 == 3] <- 0.05
tfr$sm1
tfr$sm2 <- calibrate(tfr$smb, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sm2[tfr$SM2 == 3] <- 0.05
tfr$sm2
tfr$sm3 <- calibrate(tfr$smc, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sm3[tfr$SM3 == 3] <- 0.05
tfr$sm3
tfr$sm4 <- calibrate(tfr$smd, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sm4[tfr$SM4 == 3] <- 0.05
tfr$sm4
tfr$sm5 <- calibrate(tfr$sme, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sm5[tfr$SM5 == 3] <- 0.05
tfr$sm5
tfr$sm6 <- calibrate(tfr$smf, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sm6[tfr$SM6 == 3] <- 0.05
tfr$sm6
tfr$sm7 <- calibrate(tfr$smg, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sm7[tfr$SM7 == 3] <- 0.05
tfr$sm7
tfr$sm8 <- calibrate(tfr$smh, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sm8[tfr$SM8 == 3] <- 0.05
tfr$sm8
tfr$sm9 <- calibrate(tfr$smi, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sm9[tfr$SM9 == 3] <- 0.05
tfr$sm9
tfr$sm10 <- calibrate(tfr$smj, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sm10[tfr$SM10 == 3] <- 0.05
tfr$sm10
tfr$sm11 <- calibrate(tfr$smk, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sm11[tfr$SM11 == 3] <- 0.05
tfr$sm11
tfr$sm12 <- calibrate(tfr$sml, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$sm12[tfr$SM12 == 3] <- 0.05
tfr$sm12

##For CM
## for CM
tfr$cma <- tfr$CM1
tfr$cma[tfr$CM1 == 4] <- 3
tfr$cma[tfr$CM1 == 5] <- 4
tfr$cma


tfr$cmb <- tfr$CM2
tfr$cmb[tfr$CM2 == 4] <- 3
tfr$cmb[tfr$CM2 == 5] <- 4

tfr$cmc <- tfr$CM3
tfr$cmc[tfr$CM3 == 4] <- 3
tfr$cmc[tfr$CM3 == 5] <- 4


tfr$cmd <- tfr$CM4
tfr$cmd[tfr$CM4 == 4] <- 3
tfr$cmd[tfr$CM4 == 5] <- 4

tfr$cm1 <- calibrate(tfr$cma, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$cm1[tfr$CM1 == 3] <- 0.05
tfr$cm1
tfr$cm2 <- calibrate(tfr$cmb, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$cm2[tfr$CM2 == 3] <- 0.05
tfr$cm2
tfr$cm3 <- calibrate(tfr$cmc, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$cm3[tfr$CM3 == 3] <- 0.05
tfr$cm3
tfr$cm4 <- calibrate(tfr$cmd, type = "fuzzy", thresholds = c(4, 2.5, 1), logistic = TRUE)
tfr$cm4[tfr$CM4 == 3] <- 0.05
tfr$cm4

skewcm1 <- as.numeric(tfr$cm1 > 0.5)
tabcm1 <- table(skewcm1)
prop.table(tabcm1)

skewcm2 <- as.numeric(tfr$cm2 > 0.5)
tabcm2 <- table(skewcm2)
prop.table(tabcm2)

skewcm3 <- as.numeric(tfr$cm3 > 0.5)
tabcm3 <- table(skewcm3)
prop.table(tabcm3)

skewcm4 <- as.numeric(tfr$cm4 > 0.5)
tabcm4 <- table(skewcm4)
prop.table(tabcm4)

#####Aggregation########

#build aggregated sets using the maximum rule, check for values of 0.5 and for skewedness

tfr$w <- pmax(tfr$w1, tfr$w2, tfr$w3, tfr$w4, tfr$w5, na.rm = TRUE)
tfr$w

checkw <- as.numeric(tfr$w == 0.5)
sum(checkw == 1)

wskew <- as.numeric(tfr$w > 0.5)
tabw <- table(wskew)
prop.table(tabw)

tfr$sp <- pmax(tfr$sp1, tfr$sp2, tfr$sp3, tfr$sp4, tfr$sp5, tfr$sp6, na.rm = TRUE)
tfr$sp

checksp <- as.numeric(tfr$sp == 0.5)
sum(checksp == 1)

skewsp <- as.numeric(tfr$sp > 0.5)
tabsp <- table(skewsp)
prop.table(tabsp)


tfr$tp <- pmax(tfr$tp1, tfr$tp2, tfr$tp3, tfr$tp4, tfr$tp5, tfr$tp6, na.rm = TRUE)
tfr$tp

checktp <- as.numeric(tfr$tp == 0.5)
sum(checktp == 1)

skewtp <- as.numeric(tfr$tp > 0.5)
tabtp <- table(skewtp)
prop.table(tabtp)

tfr$op <- pmax(tfr$op1, tfr$op2, tfr$op3, tfr$op4, tfr$op5, tfr$op6, na.rm = TRUE)
tfr$op

checkop <- as.numeric(tfr$op == 0.5)
sum(checkop == 1)

skewop <- as.numeric(tfr$op > 0.5)
tabop <- table(skewop)
prop.table(tabop)

tfr$sm <- pmax(tfr$sm1, tfr$sm2, tfr$sm3, tfr$sm4, tfr$sm5, tfr$sm6, 
               tfr$sm7, tfr$sm8, tfr$sm9, tfr$sm10, tfr$sm11, tfr$sm12, na.rm = TRUE)
tfr$sm

checksm <- as.numeric(tfr$sm == 0.5)
sum(checksm == 1)

skewsm <- as.numeric(tfr$sm > 0.5)
tabsm <- table(skewsm)
prop.table(tabsm)


tfr$cm <- pmax(tfr$cm1, tfr$cm2, tfr$cm3, tfr$cm4, na.rm = TRUE)
tfr$cm

checkcm <- as.numeric(tfr$cm == 0.5)
sum(checkcm == 1)

skewcm <- as.numeric(tfr$cm > 0.5)
tabcm <- table(skewcm)
prop.table(tabcm)

##check for logical AND

#build aggregated sets using the minimum rule, check for  skewedness

tfr$wm <- pmin(tfr$w1, tfr$w2, tfr$w3, tfr$w4, tfr$w5, na.rm = TRUE)

wskew <- as.numeric(tfr$wm > 0.5)
tabw <- table(wskew)
prop.table(tabw)

tfr$spm <- pmin(tfr$sp1, tfr$sp2, tfr$sp3, tfr$sp4, tfr$sp5, tfr$sp6, na.rm = TRUE)

skewsp <- as.numeric(tfr$spm > 0.5)
tabsp <- table(skewsp)
prop.table(tabsp)


tfr$tpm <- pmin(tfr$tp1, tfr$tp2, tfr$tp3, tfr$tp4, tfr$tp5, tfr$tp6, na.rm = TRUE)

skewtp <- as.numeric(tfr$tpm > 0.5)
tabtp <- table(skewtp)
prop.table(tabtp)

tfr$opm <- pmin(tfr$op1, tfr$op2, tfr$op3, tfr$op4, tfr$op5, tfr$op6, na.rm = TRUE)

skewop <- as.numeric(tfr$opm > 0.5)
tabop <- table(skewop)
prop.table(tabop)

tfr$smm <- pmin(tfr$sm1, tfr$sm2, tfr$sm3, tfr$sm4, tfr$sm5, tfr$sm6, 
               tfr$sm7, tfr$sm8, tfr$sm9, tfr$sm10, tfr$sm11, tfr$sm12, na.rm = TRUE)


skewsm <- as.numeric(tfr$smm > 0.5)
tabsm <- table(skewsm)
prop.table(tabsm)


tfr$cmm <- pmin(tfr$cm1, tfr$cm2, tfr$cm3, tfr$cm4, na.rm = TRUE)

skewcm <- as.numeric(tfr$cmm > 0.5)
tabcm <- table(skewcm)
prop.table(tabcm)

#exclude cases with set membership 0.5

tfr = tfr[tfr$w != 0.5,]
tfr = tfr[tfr$op != 0.5,] 
tfr = tfr[tfr$sp != 0.5,]
tfr = tfr[tfr$tp != 0.5,]
tfr = tfr[tfr$sm != 0.5,]
tfr = tfr[tfr$cm != 0.5,]
tfr
describe(tfr)
colnames(tfr)

# rename raw variables
tfr <- rename(tfr, W = w)
tfr <- rename(tfr, SP=  sp)
tfr <- rename(tfr, TP = tp)
tfr <- rename(tfr, OP = op)
tfr <- rename(tfr, SM = sm)
tfr <- rename(tfr, CM = cm)
colnames (tfr)


tff_d1_c1 <- subset(tfr, select = c("W", "SP", "TP", "OP", "SM", "CM"))
write.csv2(tff_d1_c1, file = "tff_d1_c1.csv")


###################ANALYSIS OF NECESSITY############

tff <- read.csv("tff_d1_c1.csv", row.names=1, header = TRUE, sep = ";",  dec = ",")
describe(tff)

#analysis of necessity for W and w  

superSubset(tff, outcome = "W", 
            conditions = "OP, SP, TP, SM, CM", 
            incl.cut = 0.9, cov.cut = 0.6)

superSubset(tff, outcome = "~W", 
            conditions = "OP, SP, TP, SM, CM", 
            incl.cut = 0.9, cov.cut = 0.6)




#calculate Z-Test for necessary conditions. I apply a RoN threshold of 0.6, so only for NCs that meet this threshold

nw <- 1-tff$W

# for NC1 OP + TP

NC1W <- compute("OP+TP", data=tff)
obsNC1W <- as.numeric(NC1W >= tff$W)

zNC1W <- round(((((sum(obsNC1W == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits = 1)
zNC1W

pzNC1W <- 2*pnorm(-abs(zNC1W))
pzNC1W

# for NC2 OP+SP

NC2W <- compute("OP+SP", data=tff)
obsNC2W <- as.numeric(NC2W >= tff$W)

zNC2W <- round(((((sum(obsNC2W == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits = 1)
zNC2W

pzNC2W <- 2*pnorm(-abs(zNC2W))
pzNC2W


# for NC3 SP+TP+SM

NC3W <- compute("SP+TP+SM", data=tff)
obsNC3W <- as.numeric(NC3W >= tff$W)

zNC3W <-round( ((((sum(obsNC3W == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zNC3W

pzNC3W <- 2*pnorm(-abs(zNC3W))
pzNC3W

# for NC4 OP+SM+CM

NC4W <- compute("OP+SM+CM ", data=tff)
obsNC4W <- as.numeric(NC4W >= tff$W)

zNC4W <-round( ((((sum(obsNC4W == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zNC4W

pzNC4W <- 2*pnorm(-abs(zNC4W))
pzNC4W



#Test for skewness of necessary conditions. It is tested whether the proportion of cases in all but one 
#specific area of the XY plot is lower than ten per cent.
##First, visual check with XY Plots; then calculate number of cases outside each area; 
#and caluclate percentage outside area; 
#there are a lot of cases on the diagonal, which are attributed to two area. This is why 
#the percentages don't add up to 100.



A1NC1W <- round(sum(as.numeric((NC1W <= tff$W) & (NC1W >= nw)))/(1004/100), digits = 1)

A2NC1W <- round(sum(as.numeric((NC1W >= tff$W)&(NC1W >= nw)))/(1004/100), digits = 1)

A3NC1W <- round(sum(as.numeric((NC1W >= tff$W)&(NC1W <= nw)))/(1004/100), digits=1)

A4NC1W <- round(sum(as.numeric((NC1W <= tff$W)&(NC1W <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 2))
plot(tff$W, NC1W, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness test for TP + OP',
     xlab='TP + OP',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1NC1W)
text(x=0.8, y=0.5, labels = A2NC1W)
text(x=0.5, y=0.2, labels = A3NC1W)
text(x=0.2, y=0.5, labels = A4NC1W)

A1NC2W <- round(sum(as.numeric((NC2W <= tff$W) & (NC2W >= nw)))/(1004/100), digits = 1)

A2NC2W <- round(sum(as.numeric((NC2W >= tff$W)&(NC2W >= nw)))/(1004/100), digits = 1)

A3NC2W <- round(sum(as.numeric((NC2W >= tff$W)&(NC2W <= nw)))/(1004/100), digits=1)

A4NC2W <- round(sum(as.numeric((NC2W <= tff$W)&(NC2W <= nw)))/(1004/100), digits = 1)

plot(tff$W, NC2W, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness test for SP + OP',
     xlab='SP + OP',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1NC2W)
text(x=0.8, y=0.5, labels = A2NC2W)
text(x=0.5, y=0.2, labels = A3NC2W)
text(x=0.2, y=0.5, labels = A4NC2W)

A1NC3W <- round(sum(as.numeric((NC3W <= tff$W) & (NC3W >= nw)))/(1004/100), digits = 1)

A2NC3W <- round(sum(as.numeric((NC3W >= tff$W)&(NC3W >= nw)))/(1004/100), digits = 1)

A3NC3W <- round(sum(as.numeric((NC3W >= tff$W)&(NC3W <= nw)))/(1004/100), digits=1)

A4NC3W <- round(sum(as.numeric((NC3W <= tff$W)&(NC3W <= nw)))/(1004/100), digits = 1)

plot(tff$W, NC3W, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness test for OP + SM + CM',
     xlab='OP + SM + CM',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1NC3W)
text(x=0.8, y=0.5, labels = A2NC3W)
text(x=0.5, y=0.2, labels = A3NC3W)
text(x=0.2, y=0.5, labels = A4NC3W)

A1NC4W <- round(sum(as.numeric((NC4W <= tff$W) & (NC4W >= nw)))/(1004/100), digits = 1)

A2NC4W <- round(sum(as.numeric((NC4W >= tff$W)&(NC4W >= nw)))/(1004/100), digits = 1)

A3NC4W <- round(sum(as.numeric((NC4W >= tff$W)&(NC4W <= nw)))/(1004/100), digits=1)

A4NC4W <- round(sum(as.numeric((NC4W <= tff$W)&(NC4W <= nw)))/(1004/100), digits = 1)
plot(tff$W, NC4W, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness test for SP + TP + SM',
     xlab='SP + TP + SM',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1NC4W)
text(x=0.8, y=0.5, labels = A2NC4W)
text(x=0.5, y=0.2, labels = A3NC4W)
text(x=0.2, y=0.5, labels = A4NC4W)


#####check accuracy of untenable assumptions####

#as the untenable assumptions for the analysis of sufficiency for W were set sp*op + tp*op, it was assumed that the two necessary conditions are connectd with the logical AND. However, with imperfect subset relations this is not always the case. 
#Check if it is the case:

necall <- compute("(TP+OP)*(SP+OP)", data=tff)

pof(necall, W, tff, relation = "nec")

# the condition has consistency necessity of .870. Hence, the untenable assumptions are slightly overly restrictive.

                      
#################ANALYSIS OF SUFFICIENCY FOR W##################
# 


tff <- read.csv("tff_d1_c1.csv", row.names=1, header = TRUE, sep = ";",  dec = ",")
describe(tff)

nw <- 1-tff$W

#analysis of sufficiency for W
# calculate untenable assumptions

nec <- deMorgan("(TP + OP)*(SP + OP)", snames = "TP, OP, SP")
nec

#### W1: frequency threshold 5, raw consistency threshold 0.924 ####

ttW1 <- truthTable(tff, outcome="W", 
                        conditions = "SP, TP, OP, SM, CM", 
                        incl.cut=0.924, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                        complete=FALSE, show.cases=FALSE)
ttW1



#####W2  frequency threshold 5, raw consistency threshold 0.922 #####


ttW2 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.922, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW2




#parsimonious solution
psW2 <- eqmcc(ttW2, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW2

#exclude contradictory assumptions & rows: tp*op + sp*op (contradict statement of necessity)

ttW2new <- ttW2
nec
rows <- findRows("tp*op + op*sp", ttW2new, remainders = FALSE)
rows

ttW2new$tt[as.character(rows), "OUT"] <- 0
ttW2new


#conservative solution
csW2 <- eqmcc(ttW2new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW2

#enhanced parsimonious solution (all simplifying assumptions)

epsW2 <- eqmcc(ttW2new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW2


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  SP*OP*SM*CM + SP*TP*SM*CM + TP*OP*SM*CM
epsolW2 <- compute("SP*OP*SM*CM + SP*TP*SM*CM + TP*OP*SM*CM", data=tff)

obsepsolW2 <- as.numeric(epsolW2 <= tff$W)
zepsolW2 <- round(((((sum(obsepsolW2 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolW2
pzepsolW2 <- 2*pnorm(-abs(zepsolW2))
pzepsolW2

#Z-Test for enhanced parsimonious solution, 0.75 (very often sufficient)
z2epsolW2 <- round(((((sum(obsepsolW2 == 1))/1004) - 0.75)-(1/(2*1004))) / (sqrt((0.75*(1-0.75))/1004)), digits=1)
z2epsolW2
p2zepsolW2 <- 2*pnorm(-abs(z2epsolW2))
p2zepsolW2


#intermediate solution with directional expectations SM -> W, CM -> W
isW2 <- eqmcc(ttW2new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW2


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M2: SP*OP*SM*CM + SP*TP*SM*CM + TP*OP*SM*CM
path1W2 <- compute("SP*OP*SM*CM ", data=tff)
path2W2 <- compute("SP*TP*SM*CM ", data=tff)
path3W2 <- compute("TP*OP*SM*CM ", data=tff)


A1path1W2 <- round(sum(as.numeric((path1W2 <= tff$W) & (path1W2 >= nw)))/(1004/100), digits = 1)

A2path1W2 <- round(sum(as.numeric((path1W2 >= tff$W)&(path1W2 >= nw)))/(1004/100), digits = 1)

A3path1W2 <- round(sum(as.numeric((path1W2 >= tff$W)&(path1W2 <= nw)))/(1004/100), digits=1)

A4path1W2 <- round(sum(as.numeric((path1W2 <= tff$W)&(path1W2 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W2)
text(x=0.8, y=0.5, labels = A2path1W2)
text(x=0.5, y=0.2, labels = A3path1W2)
text(x=0.2, y=0.5, labels = A4path1W2)



A1path2W2 <- round(sum(as.numeric((path2W2 <= tff$W) & (path2W2 >= nw)))/(1004/100), digits = 1)

A2path2W2 <- round(sum(as.numeric((path2W2 >= tff$W)&(path2W2 >= nw)))/(1004/100), digits = 1)

A3path2W2 <- round(sum(as.numeric((path2W2 >= tff$W)&(path2W2 <= nw)))/(1004/100), digits=1)

A4path2W2 <- round(sum(as.numeric((path2W2 <= tff$W)&(path2W2 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path2W2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab='Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W2)
text(x=0.8, y=0.5, labels = A2path2W2)
text(x=0.5, y=0.2, labels = A3path2W2)
text(x=0.2, y=0.5, labels = A4path2W2)

A1path3W2 <- round(sum(as.numeric((path3W2 <= tff$W) & (path3W2 >= nw)))/(1004/100), digits = 1)

A2path3W2 <- round(sum(as.numeric((path3W2 >= tff$W)&(path3W2 >= nw)))/(1004/100), digits = 1)

A3path3W2 <- round(sum(as.numeric((path3W2 >= tff$W)&(path3W2 <= nw)))/(1004/100), digits=1)

A4path3W2 <- round(sum(as.numeric((path3W2 <= tff$W)&(path3W2 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path3W2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness path 3',
     xlab='Path 3',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3W2)
text(x=0.8, y=0.5, labels = A2path3W2)
text(x=0.5, y=0.2, labels = A3path3W2)
text(x=0.2, y=0.5, labels = A4path3W2)



A1epsolW2 <- round(sum(as.numeric((epsolW2 <= tff$W) & (epsolW2 >= nw)))/(1004/100), digits = 1)

A2epsolW2 <- round(sum(as.numeric((epsolW2 >= tff$W)&(epsolW2 >= nw)))/(1004/100), digits = 1)

A3epsolW2 <- round(sum(as.numeric((epsolW2 >= tff$W)&(epsolW2 <= nw)))/(1004/100), digits=1)

A4epsolW2 <- round(sum(as.numeric((epsolW2 <= tff$W)&(epsolW2 <= nw)))/(1004/100), digits = 1)

isW2
plot(tff$W, epsolW2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW2)
text(x=0.8, y=0.5, labels = A2epsolW2)
text(x=0.5, y=0.2, labels = A3epsolW2)
text(x=0.2, y=0.5, labels = A4epsolW2)


####W3: frequency threshold 5, raw consistency threshold 0.910#######

ttW3 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.910, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW3


#parsimonious solution
psW3 <- eqmcc(ttW3, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW3

#exclude contradictory assumptions & rows: tp*op + sp*op (contradict statement of necessity)

ttW3new <- ttW3
nec
rows <- findRows("tp*op + op*sp", ttW3new, remainders = FALSE)
rows

ttW3new$tt[as.character(rows), "OUT"] <- 0
ttW3new


#conservative solution
csW3 <- eqmcc(ttW3new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW3


#enhanced parsimonious solution (all simplifying assumptions)

epsW3 <- eqmcc(ttW3new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW3

#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term OP*SM*CM + SP*TP*SM*CM
epsolW3 <- compute("OP*SM*CM + SP*TP*SM*CM", data=tff)

obsepsolW3 <- as.numeric(epsolW3 <= tff$W)
zepsolW3 <- round(((((sum(obsepsolW3 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolW3
pzepsolW3 <- 2*pnorm(-abs(zepsolW3))
pzepsolW3

#Z-Test for enhanced parsimonious solution, 0.75 (very often sufficient)
z2epsolW3 <- round(((((sum(obsepsolW3 == 1))/1004) - 0.75)-(1/(2*1004))) / (sqrt((0.75*(1-0.75))/1004)), digits=1)
z2epsolW3
p2zepsolW3 <- 2*pnorm(-abs(z2epsolW3))
p2zepsolW3

#intermediate solution with directional expectations SM -> W, CM -> W
isW3 <- eqmcc(ttW3new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW3

# identify simplifying assumptions and easy counterfactuals

SAW3 <- psW3$SA
SAW3

ECW3 <- isW3$i.sol$C1P1$EC
ECW3

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: OP*SM*CM + SP*TP*SM*CM

path1W3 <- compute("OP*SM*CM", data=tff)
path2W3 <- compute("SP*TP*SM*CM ", data=tff)


A1path1W3 <- round(sum(as.numeric((path1W3 <= tff$W) & (path1W3 >= nw)))/(1004/100), digits = 1)

A2path1W3 <- round(sum(as.numeric((path1W3 >= tff$W)&(path1W3 >= nw)))/(1004/100), digits = 1)

A3path1W3 <- round(sum(as.numeric((path1W3 >= tff$W)&(path1W3 <= nw)))/(1004/100), digits=1)

A4path1W3 <- round(sum(as.numeric((path1W3 <= tff$W)&(path1W3 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W3)
text(x=0.8, y=0.5, labels = A2path1W3)
text(x=0.5, y=0.2, labels = A3path1W3)
text(x=0.2, y=0.5, labels = A4path1W3)



A1path2W3 <- round(sum(as.numeric((path2W3 <= tff$W) & (path2W3 >= nw)))/(1004/100), digits = 1)

A2path2W3 <- round(sum(as.numeric((path2W3 >= tff$W)&(path2W3 >= nw)))/(1004/100), digits = 1)

A3path2W3 <- round(sum(as.numeric((path2W3 >= tff$W)&(path2W3 <= nw)))/(1004/100), digits=1)

A4path2W3 <- round(sum(as.numeric((path2W3 <= tff$W)&(path2W3 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path2W3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab='Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W3)
text(x=0.8, y=0.5, labels = A2path2W3)
text(x=0.5, y=0.2, labels = A3path2W3)
text(x=0.2, y=0.5, labels = A4path2W3)



A1epsolW3 <- round(sum(as.numeric((epsolW3 <= tff$W) & (epsolW3 >= nw)))/(1004/100), digits = 1)

A2epsolW3 <- round(sum(as.numeric((epsolW3 >= tff$W)&(epsolW3 >= nw)))/(1004/100), digits = 1)

A3epsolW3 <- round(sum(as.numeric((epsolW3 >= tff$W)&(epsolW3 <= nw)))/(1004/100), digits=1)

A4epsolW3 <- round(sum(as.numeric((epsolW3 <= tff$W)&(epsolW3 <= nw)))/(1004/100), digits = 1)


plot(tff$W, epsolW3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW3)
text(x=0.8, y=0.5, labels = A2epsolW3)
text(x=0.5, y=0.2, labels = A3epsolW3)
text(x=0.2, y=0.5, labels = A4epsolW3)


#Z-Test for intermediate solution, 0.75 (very often sufficient)
# Calculate membership in solution term  OP*SM*CM + SP*TP*SM*CM
isW3
isolW3 <- compute("OP*SM*CM + SP*TP*SM*CM", data=tff)

obsisolW3 <- as.numeric(isolW3 <= tff$W)
zisolW3 <- round(((((sum(obsisolW3 == 1))/1004) - 0.75)-(1/(2*1004))) / (sqrt((0.75*(1-0.75))/1004)), digits=1)
zisolW3
pzisolW3 <- 2*pnorm(-abs(zisolW3))
pzisolW3


# for path 1
# Calculate membership in path  OP*SM*CM
path1isolW3 <- compute("OP*SM*CM ", data=tff)

obspath1isolW3 <- as.numeric(path1isolW3 <= tff$W)
zpath1isolW3 <- round(((((sum(obspath1isolW3 == 1))/1004) - 0.75)-(1/(2*1004))) / (sqrt((0.75*(1-0.75))/1004)), digits=1)
zpath1isolW3
pzpath1isolW3 <- 2*pnorm(-abs(zpath1isolW3))
pzpath1isolW3

#for path 2
# Calculate membership in path  SP*TP*SM*CM
path2isolW3 <- compute("SP*TP*SM*CM ", data=tff)

obspath2isolW3 <- as.numeric(path2isolW3 <= tff$W)
zpath2isolW3 <- round(((((sum(obspath2isolW3 == 1))/1004) - 0.75)-(1/(2*1004))) / (sqrt((0.75*(1-0.75))/1004)), digits=1)
zpath2isolW3
pzpath2isolW3 <- 2*pnorm(-abs(zpath2isolW3))
pzpath2isolW3


####W4 frequency threshold 10, raw consistency threshold 0.906####


ttW4 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.906, n.cut=10, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW4


colnames(tff)

#parsimonious solution
psW4 <- eqmcc(ttW4, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW4

#exclude contradictory assumptions & rows: tp*op + sp*op (contradict statement of necessity)

ttW4new <- ttW4
nec
rows <- findRows("tp*op + op*sp", ttW4new, remainders = FALSE)
rows

ttW4new$tt[as.character(rows), "OUT"] <- 0
ttW4new


#conservative solution
csW4 <- eqmcc(ttW4new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW4

#enhanced parsimonious solution (all simplifying assumptions)

epsW4 <- eqmcc(ttW4new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW4


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term OP*SM*CM + sp*TP*OP*SM
epsolW4 <- compute("OP*SM*CM + sp*TP*OP*SM", data=tff)

obsepsolW4 <- as.numeric(epsolW4 <= tff$W)
zepsolW4 <- round(((((sum(obsepsolW4 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolW4
pzepsolW4 <- 2*pnorm(-abs(zepsolW4))
pzepsolW4

#Z-Test for enhanced parsimonious solution, 0.75 (very often sufficient)
z2epsolW4 <- round(((((sum(obsepsolW4 == 1))/1004) - 0.75)-(1/(2*1004))) / (sqrt((0.75*(1-0.75))/1004)), digits=1)
z2epsolW4
p2zepsolW4 <- 2*pnorm(-abs(z2epsolW4))
p2zepsolW4

#intermediate solution with directional expectations SM -> W, CM -> W
isW4 <- eqmcc(ttW4new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW4

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: OP*SM*CM + sp*TP*OP*SM

path1W4 <- compute("OP*SM*CM ", data=tff)
path2W4 <- compute("sp*TP*OP*SM ", data=tff)


A1path1W4 <- round(sum(as.numeric((path1W4 <= tff$W) & (path1W4 >= nw)))/(1004/100), digits = 1)

A2path1W4 <- round(sum(as.numeric((path1W4 >= tff$W)&(path1W4 >= nw)))/(1004/100), digits = 1)

A3path1W4 <- round(sum(as.numeric((path1W4 >= tff$W)&(path1W4 <= nw)))/(1004/100), digits=1)

A4path1W4 <- round(sum(as.numeric((path1W4 <= tff$W)&(path1W4 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W4)
text(x=0.8, y=0.5, labels = A2path1W4)
text(x=0.5, y=0.2, labels = A3path1W4)
text(x=0.2, y=0.5, labels = A4path1W4)



A1path2W4 <- round(sum(as.numeric((path2W4 <= tff$W) & (path2W4 >= nw)))/(1004/100), digits = 1)

A2path2W4 <- round(sum(as.numeric((path2W4 >= tff$W)&(path2W4 >= nw)))/(1004/100), digits = 1)

A3path2W4 <- round(sum(as.numeric((path2W4 >= tff$W)&(path2W4 <= nw)))/(1004/100), digits=1)

A4path2W4 <- round(sum(as.numeric((path2W4 <= tff$W)&(path2W4 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path2W4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab='Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W4)
text(x=0.8, y=0.5, labels = A2path2W4)
text(x=0.5, y=0.2, labels = A3path2W4)
text(x=0.2, y=0.5, labels = A4path2W4)


A1epsolW4 <- round(sum(as.numeric((epsolW4 <= tff$W) & (epsolW4 >= nw)))/(1004/100), digits = 1)

A2epsolW4 <- round(sum(as.numeric((epsolW4 >= tff$W)&(epsolW4 >= nw)))/(1004/100), digits = 1)

A3epsolW4 <- round(sum(as.numeric((epsolW4 >= tff$W)&(epsolW4 <= nw)))/(1004/100), digits=1)

A4epsolW4 <- round(sum(as.numeric((epsolW4 <= tff$W)&(epsolW4 <= nw)))/(1004/100), digits = 1)


plot(tff$W, epsolW4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW4)
text(x=0.8, y=0.5, labels = A2epsolW4)
text(x=0.5, y=0.2, labels = A3epsolW4)
text(x=0.2, y=0.5, labels = A4epsolW4)

###### W5: frequency threshold 10, raw consistency threshold 0.900######

ttW5 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.900, n.cut=10, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW5



#parsimonious solution
psW5 <- eqmcc(ttW5, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW5

#exclude contradictory assumptions & rows: tp*op + sp*op (contradict statement of necessity)

ttW5new <- ttW5
nec
rows <- findRows("tp*op + op*sp", ttW5new, remainders = FALSE)
rows

ttW5new$tt[as.character(rows), "OUT"] <- 0
ttW5new


#conservative solution
csW5 <- eqmcc(ttW5new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW5

#enhanced parsimonious solution (all simplifying assumptions)

epsW5 <- eqmcc(ttW5new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW5

#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term OP*SM*CM + sp*TP*OP*SM
epsolW5 <- compute("OP*SM*CM + sp*TP*OP*SM", data=tff)

obsepsolW5 <- as.numeric(epsolW5 <= tff$W)
zepsolW5 <- round(((((sum(obsepsolW5 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolW5
pzepsolW5 <- 2*pnorm(-abs(zepsolW5))
pzepsolW5

#Z-Test for enhanced parsimonious solution, 0.75 (very often sufficient)
z2epsolW5 <- round(((((sum(obsepsolW5 == 1))/1004) - 0.75)-(1/(2*1004))) / (sqrt((0.75*(1-0.75))/1004)), digits=1)
z2epsolW5
p2zepsolW5 <- 2*pnorm(-abs(z2epsolW5))
p2zepsolW5

#intermediate solution with directional expectations SM -> W, CM -> W
isW5 <- eqmcc(ttW5new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW5



#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: OP*CM + sp*TP*OP*SM

path1W5 <- compute("OP*SM*CM ", data=tff)
path2W5 <- compute("sp*TP*OP*SM ", data=tff)

A1path1W5 <- round(sum(as.numeric((path1W5 <= tff$W) & (path1W5 >= nw)))/(1004/100), digits = 1)

A2path1W5 <- round(sum(as.numeric((path1W5 >= tff$W)&(path1W5 >= nw)))/(1004/100), digits = 1)

A3path1W5 <- round(sum(as.numeric((path1W5 >= tff$W)&(path1W5 <= nw)))/(1004/100), digits=1)

A4path1W5 <- round(sum(as.numeric((path1W5 <= tff$W)&(path1W5 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W5)
text(x=0.8, y=0.5, labels = A2path1W5)
text(x=0.5, y=0.2, labels = A3path1W5)
text(x=0.2, y=0.5, labels = A4path1W5)



A1path2W5 <- round(sum(as.numeric((path2W5 <= tff$W) & (path2W5 >= nw)))/(1004/100), digits = 1)

A2path2W5 <- round(sum(as.numeric((path2W5 >= tff$W)&(path2W5 >= nw)))/(1004/100), digits = 1)

A3path2W5 <- round(sum(as.numeric((path2W5 >= tff$W)&(path2W5 <= nw)))/(1004/100), digits=1)

A4path2W5 <- round(sum(as.numeric((path2W5 <= tff$W)&(path2W5 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path2W5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab='Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W5)
text(x=0.8, y=0.5, labels = A2path2W5)
text(x=0.5, y=0.2, labels = A3path2W5)
text(x=0.2, y=0.5, labels = A4path2W5)

A1epsolW5 <- round(sum(as.numeric((epsolW5 <= tff$W) & (epsolW5 >= nw)))/(1004/100), digits = 1)

A2epsolW5 <- round(sum(as.numeric((epsolW5 >= tff$W)&(epsolW5 >= nw)))/(1004/100), digits = 1)

A3epsolW5 <- round(sum(as.numeric((epsolW5 >= tff$W)&(epsolW5 <= nw)))/(1004/100), digits=1)

A4epsolW5 <- round(sum(as.numeric((epsolW5 <= tff$W)&(epsolW5 <= nw)))/(1004/100), digits = 1)

isW5
plot(tff$W, epsolW5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW5)
text(x=0.8, y=0.5, labels = A2epsolW5)
text(x=0.5, y=0.2, labels = A3epsolW5)
text(x=0.2, y=0.5, labels = A4epsolW5)



####W6 frequency threshold 10, raw consistency threshold 0.898####
ttW6 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.898, n.cut=10, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW6


colnames(tff)

#parsimonious solution
psW6 <- eqmcc(ttW6, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW6

#exclude contradictory assumptions & rows: tp*op + sp*op (contradict statement of necessity)

ttW6new <- ttW6
nec
rows <- findRows("tp*op + op*sp", ttW6new, remainders = FALSE)
rows

ttW6new$tt[as.character(rows), "OUT"] <- 0
ttW6new


#conservative solution
csW6 <- eqmcc(ttW6new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW6

#enhanced parsimonious solution (all simplifying assumptions)

epsW6 <- eqmcc(ttW6new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW6


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term OP*SM*CM + sp*TP*OP*SM + tp*OP*CM
epsolW6 <- compute("OP*SM*CM + sp*TP*OP*SM + tp*OP*CM", data=tff)

obsepsolW6 <- as.numeric(epsolW6 <= tff$W)
zepsolW6 <- round(((((sum(obsepsolW6 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolW6
pzepsolW6 <- 2*pnorm(-abs(zepsolW6))
pzepsolW6

#Z-Test for enhanced parsimonious solution, 0.75 (very often sufficient)
z2epsolW6 <- round(((((sum(obsepsolW6 == 1))/1004) - 0.75)-(1/(2*1004))) / (sqrt((0.75*(1-0.75))/1004)), digits=1)
z2epsolW6
p2zepsolW6 <- 2*pnorm(-abs(z2epsolW6))
p2zepsolW6

#intermediate solution with directional expectations SM -> W, CM -> W
isW6 <- eqmcc(ttW6new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW6

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M2: OP*SM*CM + sp*TP*OP*SM + tp*OP*CM

path1W6 <- compute("OP*SM*CM ", data=tff)
path2W6 <- compute("sp*TP*OP*SM ", data=tff)
path3W6 <- compute("tp*OP*CM ", data=tff)

A1path1W6 <- round(sum(as.numeric((path1W6 <= tff$W) & (path1W6 >= nw)))/(1004/100), digits = 1)

A2path1W6 <- round(sum(as.numeric((path1W6 >= tff$W)&(path1W6 >= nw)))/(1004/100), digits = 1)

A3path1W6 <- round(sum(as.numeric((path1W6 >= tff$W)&(path1W6 <= nw)))/(1004/100), digits=1)

A4path1W6 <- round(sum(as.numeric((path1W6 <= tff$W)&(path1W6 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W6)
text(x=0.8, y=0.5, labels = A2path1W6)
text(x=0.5, y=0.2, labels = A3path1W6)
text(x=0.2, y=0.5, labels = A4path1W6)



A1path2W6 <- round(sum(as.numeric((path2W6 <= tff$W) & (path2W6 >= nw)))/(1004/100), digits = 1)

A2path2W6 <- round(sum(as.numeric((path2W6 >= tff$W)&(path2W6 >= nw)))/(1004/100), digits = 1)

A3path2W6 <- round(sum(as.numeric((path2W6 >= tff$W)&(path2W6 <= nw)))/(1004/100), digits=1)

A4path2W6 <- round(sum(as.numeric((path2W6 <= tff$W)&(path2W6 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path2W6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab='Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W6)
text(x=0.8, y=0.5, labels = A2path2W6)
text(x=0.5, y=0.2, labels = A3path2W6)
text(x=0.2, y=0.5, labels = A4path2W6)


A1path3W6 <- round(sum(as.numeric((path3W6 <= tff$W) & (path3W6 >= nw)))/(1004/100), digits = 1)

A2path3W6 <- round(sum(as.numeric((path3W6 >= tff$W)&(path3W6 >= nw)))/(1004/100), digits = 1)

A3path3W6 <- round(sum(as.numeric((path3W6 >= tff$W)&(path3W6 <= nw)))/(1004/100), digits=1)

A4path3W6 <- round(sum(as.numeric((path3W6 <= tff$W)&(path3W6 <= nw)))/(1004/100), digits = 1)

plot(tff$W, path3W6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 3',
     xlab='Path 3',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3W6)
text(x=0.8, y=0.5, labels = A2path3W6)
text(x=0.5, y=0.2, labels = A3path3W6)
text(x=0.2, y=0.5, labels = A4path3W6)

A1epsolW6 <- round(sum(as.numeric((epsolW6 <= tff$W) & (epsolW6 >= nw)))/(1004/100), digits = 1)

A2epsolW6 <- round(sum(as.numeric((epsolW6 >= tff$W)&(epsolW6 >= nw)))/(1004/100), digits = 1)

A3epsolW6 <- round(sum(as.numeric((epsolW6 >= tff$W)&(epsolW6 <= nw)))/(1004/100), digits=1)

A4epsolW6 <- round(sum(as.numeric((epsolW6 <= tff$W)&(epsolW6 <= nw)))/(1004/100), digits = 1)

isW6
plot(tff$W, epsolW6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW6)
text(x=0.8, y=0.5, labels = A2epsolW6)
text(x=0.5, y=0.2, labels = A3epsolW6)
text(x=0.2, y=0.5, labels = A4epsolW6)


#####W7 frequency threshold 50, raw consistency threshold 0.824####

ttW7 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.824, n.cut=50, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW7


colnames(tff)

#parsimonious solution
psW7 <- eqmcc(ttW7, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW7

#exclude contradictory assumptions & rows: tp*op + sp*op (contradict statement of necessity)

ttW7new <- ttW7
nec
rows <- findRows("tp*op + op*sp", ttW7new, remainders = FALSE)
rows

ttW7new$tt[as.character(rows), "OUT"] <- 0
ttW7new


#conservative solution
csW7 <- eqmcc(ttW7new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW7

#enhanced parsimonious solution (all simplifying assumptions)

epsW7 <- eqmcc(ttW7new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW7


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term SP*TP*SM
epsolW7 <- compute("SP*TP*SM", data=tff)

obsepsolW7 <- as.numeric(epsolW7 <= tff$W)
zepsolW7 <- round(((((sum(obsepsolW7 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolW7
pzepsolW7 <- 2*pnorm(-abs(zepsolW7))
pzepsolW7

#Z-Test for enhanced parsimonious solution, 0.75 (very often sufficient)
z2epsolW7 <- round(((((sum(obsepsolW7 == 1))/1004) - 0.75)-(1/(2*1004))) / (sqrt((0.75*(1-0.75))/1004)), digits=1)
z2epsolW7
p2zepsolW7 <- 2*pnorm(-abs(z2epsolW7))
p2zepsolW7

#intermediate solution with directional expectations SM -> W, CM -> W
isW7 <- eqmcc(ttW7new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW7

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M2: SP*TP*SM

path1W7 <- compute("SP*TP*SM", data=tff)

A1path1W7 <- round(sum(as.numeric((path1W7 <= tff$W) & (path1W7 >= nw)))/(1004/100), digits = 1)

A2path1W7 <- round(sum(as.numeric((path1W7 >= tff$W)&(path1W7 >= nw)))/(1004/100), digits = 1)

A3path1W7 <- round(sum(as.numeric((path1W7 >= tff$W)&(path1W7 <= nw)))/(1004/100), digits=1)

A4path1W7 <- round(sum(as.numeric((path1W7 <= tff$W)&(path1W7 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W7, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W7)
text(x=0.8, y=0.5, labels = A2path1W7)
text(x=0.5, y=0.2, labels = A3path1W7)
text(x=0.2, y=0.5, labels = A4path1W7)


A1epsolW7 <- round(sum(as.numeric((epsolW7 <= tff$W) & (epsolW7 >= nw)))/(1004/100), digits = 1)

A2epsolW7 <- round(sum(as.numeric((epsolW7 >= tff$W)&(epsolW7 >= nw)))/(1004/100), digits = 1)

A3epsolW7 <- round(sum(as.numeric((epsolW7 >= tff$W)&(epsolW7 <= nw)))/(1004/100), digits=1)

A4epsolW7 <- round(sum(as.numeric((epsolW7 <= tff$W)&(epsolW7 <= nw)))/(1004/100), digits = 1)

isW7
plot(tff$W, epsolW7, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW7)
text(x=0.8, y=0.5, labels = A2epsolW7)
text(x=0.5, y=0.2, labels = A3epsolW7)
text(x=0.2, y=0.5, labels = A4epsolW7)

#####W8 frequency threshold 50, raw consistency threshold 0.802####

ttW8 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.802, n.cut=50, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW8


colnames(tff)

#parsimonious solution
psW8 <- eqmcc(ttW8, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW8

#exclude contradictory assumptions & rows: tp*op + sp*op (contradict statement of necessity)

ttW8new <- ttW8
nec
rows <- findRows("tp*op + op*sp", ttW8new, remainders = FALSE)
rows

ttW8new$tt[as.character(rows), "OUT"] <- 0
ttW8new


#conservative solution
csW8 <- eqmcc(ttW8new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW8

#enhanced parsimonious solution (all simplifying assumptions)

epsW8 <- eqmcc(ttW8new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW8


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term SP*OP
epsolW8 <- compute("SP*OP", data=tff)

obsepsolW8 <- as.numeric(epsolW8 <= tff$W)
zepsolW8 <- round(((((sum(obsepsolW8 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolW8
pzepsolW8 <- 2*pnorm(-abs(zepsolW8))
pzepsolW8

#Z-Test for enhanced parsimonious solution, 0.75 (very often sufficient)
z2epsolW8 <- round(((((sum(obsepsolW8 == 1))/1004) - 0.75)-(1/(2*1004))) / (sqrt((0.75*(1-0.75))/1004)), digits=1)
z2epsolW8
p2zepsolW8 <- 2*pnorm(-abs(z2epsolW8))
p2zepsolW8


#intermediate solution with directional expectations SM -> W, CM -> W
isW8 <- eqmcc(ttW8new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW8

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: SP*OP

path1W8 <- compute("SP*OP ", data=tff)

A1path1W8 <- round(sum(as.numeric((path1W8 <= tff$W) & (path1W8 >= nw)))/(1004/100), digits = 1)

A2path1W8 <- round(sum(as.numeric((path1W8 >= tff$W)&(path1W8 >= nw)))/(1004/100), digits = 1)

A3path1W8 <- round(sum(as.numeric((path1W8 >= tff$W)&(path1W8 <= nw)))/(1004/100), digits=1)

A4path1W8 <- round(sum(as.numeric((path1W8 <= tff$W)&(path1W8 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W8, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W8)
text(x=0.8, y=0.5, labels = A2path1W8)
text(x=0.5, y=0.2, labels = A3path1W8)
text(x=0.2, y=0.5, labels = A4path1W8)


A1epsolW8 <- round(sum(as.numeric((epsolW8 <= tff$W) & (epsolW8 >= nw)))/(1004/100), digits = 1)

A2epsolW8 <- round(sum(as.numeric((epsolW8 >= tff$W)&(epsolW8 >= nw)))/(1004/100), digits = 1)

A3epsolW8 <- round(sum(as.numeric((epsolW8 >= tff$W)&(epsolW8 <= nw)))/(1004/100), digits=1)

A4epsolW8 <- round(sum(as.numeric((epsolW8 <= tff$W)&(epsolW8 <= nw)))/(1004/100), digits = 1)

isW8
plot(tff$W, epsolW8, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW8)
text(x=0.8, y=0.5, labels = A2epsolW8)
text(x=0.5, y=0.2, labels = A3epsolW8)
text(x=0.2, y=0.5, labels = A4epsolW8)

#####W9 frequency threshold 50, raw consistency threshold 0.799####

ttW9 <- truthTable(tff, outcome="W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.799, n.cut=50, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttW9


#parsimonious solution
psW9 <- eqmcc(ttW9, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psW9

#exclude contradictory assumptions & rows: tp*op + sp*op (contradict statement of necessity)

ttW9new <- ttW9
nec
rows <- findRows("tp*op + op*sp", ttW9new, remainders = FALSE)
rows

ttW9new$tt[as.character(rows), "OUT"] <- 0
ttW9new


#conservative solution
csW9 <- eqmcc(ttW9new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csW9

#enhanced parsimonious solution (all simplifying assumptions)

epsW9 <- eqmcc(ttW9new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsW9

#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term SP*OP + TP*OP
epsolW9 <- compute("SP*OP + TP*OP", data=tff)

obsepsolW9 <- as.numeric(epsolW9 <= tff$W)
zepsolW9 <- round(((((sum(obsepsolW9 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolW9
pzepsolW9 <- 2*pnorm(-abs(zepsolW9))
pzepsolW9

#Z-Test for enhanced parsimonious solution, 0.75 (very often sufficient)
z2epsolW9 <- round(((((sum(obsepsolW9 == 1))/1004) - 0.75)-(1/(2*1004))) / (sqrt((0.75*(1-0.75))/1004)), digits=1)
z2epsolW9
p2zepsolW9 <- 2*pnorm(-abs(z2epsolW9))
p2zepsolW9

#intermediate solution with directional expectations SM -> W, CM -> W
isW9 <- eqmcc(ttW9new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "-, -, -, 1, 1")

isW9

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: SP*OP + TP*OP

path1W9 <- compute("SP*OP ", data=tff)

A1path1W9 <- round(sum(as.numeric((path1W9 <= tff$W) & (path1W9 >= nw)))/(1004/100), digits = 1)

A2path1W9 <- round(sum(as.numeric((path1W9 >= tff$W)&(path1W9 >= nw)))/(1004/100), digits = 1)

A3path1W9 <- round(sum(as.numeric((path1W9 >= tff$W)&(path1W9 <= nw)))/(1004/100), digits=1)

A4path1W9 <- round(sum(as.numeric((path1W9 <= tff$W)&(path1W9 <= nw)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(tff$W, path1W9, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1W9)
text(x=0.8, y=0.5, labels = A2path1W9)
text(x=0.5, y=0.2, labels = A3path1W9)
text(x=0.2, y=0.5, labels = A4path1W9)

path2W9 <- compute("TP*OP ", data=tff)

A1path2W9 <- round(sum(as.numeric((path2W9 <= tff$W) & (path2W9 >= nw)))/(1004/100), digits = 1)

A2path2W9 <- round(sum(as.numeric((path2W9 >= tff$W)&(path2W9 >= nw)))/(1004/100), digits = 1)

A3path2W9 <- round(sum(as.numeric((path2W9 >= tff$W)&(path2W9 <= nw)))/(1004/100), digits=1)

A4path2W9 <- round(sum(as.numeric((path2W9 <= tff$W)&(path2W9 <= nw)))/(1004/100), digits = 1)


plot(tff$W, path2W9, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2W9)
text(x=0.8, y=0.5, labels = A2path2W9)
text(x=0.5, y=0.2, labels = A3path2W9)
text(x=0.2, y=0.5, labels = A4path2W9)



A1epsolW9 <- round(sum(as.numeric((epsolW9 <= tff$W) & (epsolW9 >= nw)))/(1004/100), digits = 1)

A2epsolW9 <- round(sum(as.numeric((epsolW9 >= tff$W)&(epsolW9 >= nw)))/(1004/100), digits = 1)

A3epsolW9 <- round(sum(as.numeric((epsolW9 >= tff$W)&(epsolW9 <= nw)))/(1004/100), digits=1)

A4epsolW9 <- round(sum(as.numeric((epsolW9 <= tff$W)&(epsolW9 <= nw)))/(1004/100), digits = 1)

isW9
plot(tff$W, epsolW9, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main=' Skewness eps',
     xlab='Solution term',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolW9)
text(x=0.8, y=0.5, labels = A2epsolW9)
text(x=0.5, y=0.2, labels = A3epsolW9)
text(x=0.2, y=0.5, labels = A4epsolW9)




#######################ANALYSIS OF SUFFICIENCY FOR w ####################


tff <- read.csv("tff_d1_c1.csv", row.names=1, header = TRUE, sep = ";",  dec = ",")
describe(tff)

nw <- 1-tff$W

# unenable assumptions: 

suf <- "OP*SM*CM + SP*TP*SM*CM"
suf



#####w1  frequency threshold 5, raw consistency threshold 0.972 #####


ttw1 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.972, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw1


#parsimonious solution
psw1 <- eqmcc(ttw1, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw1

#exclude contradictory assumptions & rows: OP*SM*CM + SP*TP*SM*CM  (contradict statement of sufficiency)

ttw1new <- ttw1
suf
rows <- findRows("OP*SM*CM + SP*TP*SM*CM", ttw1new, remainders = FALSE)
rows

ttw1new$tt[as.character(rows), "OUT"] <- 0
ttw1new


#conservative solution
csw1 <- eqmcc(ttw1new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw1

#enhanced parsimonious solution (all simplifying assumptions)

epsw1 <- eqmcc(ttw1new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw1


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  op*sm*CM + sp*op*CM + tp*sm*CM + SP*op*SM*cm
epsolw1 <- compute("op*sm*CM + sp*op*CM + tp*sm*CM + SP*op*SM*cm ", data=tff)

obsepsolw1 <- as.numeric(epsolw1 <= nw)
zepsolw1 <- round(((((sum(obsepsolw1 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolw1
pzepsolw1 <- 2*pnorm(-abs(zepsolw1))
pzepsolw1


#intermediate solution with Directional expectations: sp->w, tp->w, op->w.

isw1 <- eqmcc(ttw1new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw1


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: op*sm*CM + sp*op*CM + tp*sm*CM + SP*op*SM*cm
path1w1 <- compute("op*sm*CM ", data=tff)
path2w1 <- compute("sp*op*CM ", data=tff)
path3w1 <- compute("tp*sm*CM ", data=tff)
path4w1 <- compute("SP*op*SM*cm ", data=tff)

A1path1w1 <- round(sum(as.numeric((path1w1 <= nw) & (path1w1 >= tff$W)))/(1004/100), digits = 1)

A2path1w1 <- round(sum(as.numeric((path1w1 >= nw)&(path1w1 >= tff$W)))/(1004/100), digits = 1)

A3path1w1 <- round(sum(as.numeric((path1w1 >= nw)&(path1w1 <= tff$W)))/(1004/100), digits=1)

A4path1w1 <- round(sum(as.numeric((path1w1 <= nw)&(path1w1 <= tff$W)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w1)
text(x=0.8, y=0.5, labels = A2path1w1)
text(x=0.5, y=0.2, labels = A3path1w1)
text(x=0.2, y=0.5, labels = A4path1w1)



A1path2w1 <- round(sum(as.numeric((path2w1 <= nw) & (path2w1 >= tff$W)))/(1004/100), digits = 1)

A2path2w1 <- round(sum(as.numeric((path2w1 >= nw)&(path2w1 >= tff$W)))/(1004/100), digits = 1)

A3path2w1 <- round(sum(as.numeric((path2w1 >= nw)&(path2w1 <= tff$W)))/(1004/100), digits=1)

A4path2w1 <- round(sum(as.numeric((path2w1 <= nw)&(path2w1 <= tff$W)))/(1004/100), digits = 1)


plot(nw, path2w1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w1)
text(x=0.8, y=0.5, labels = A2path2w1)
text(x=0.5, y=0.2, labels = A3path2w1)
text(x=0.2, y=0.5, labels = A4path2w1)


A1path3w1 <- round(sum(as.numeric((path3w1 <= nw) & (path3w1 >= tff$W)))/(1004/100), digits = 1)

A2path3w1 <- round(sum(as.numeric((path3w1 >= nw)&(path3w1 >= tff$W)))/(1004/100), digits = 1)

A3path3w1 <- round(sum(as.numeric((path3w1 >= nw)&(path3w1 <= tff$W)))/(1004/100), digits=1)

A4path3w1 <- round(sum(as.numeric((path3w1 <= nw)&(path3w1 <= tff$W)))/(1004/100), digits = 1)


plot(nw, path3w1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 3',
     xlab= 'Path 3',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3w1)
text(x=0.8, y=0.5, labels = A2path3w1)
text(x=0.5, y=0.2, labels = A3path3w1)
text(x=0.2, y=0.5, labels = A4path3w1)


A1path4w1 <- round(sum(as.numeric((path4w1 <= nw) & (path4w1 >= tff$W)))/(1004/100), digits = 1)

A2path4w1 <- round(sum(as.numeric((path4w1 >= nw)&(path4w1 >= tff$W)))/(1004/100), digits = 1)

A3path4w1 <- round(sum(as.numeric((path4w1 >= nw)&(path4w1 <= tff$W)))/(1004/100), digits=1)

A4path4w1 <- round(sum(as.numeric((path4w1 <= nw)&(path4w1 <= tff$W)))/(1004/100), digits = 1)


plot(nw, path4w1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 4',
     xlab= 'Path 4',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path4w1)
text(x=0.8, y=0.5, labels = A2path4w1)
text(x=0.5, y=0.2, labels = A3path4w1)
text(x=0.2, y=0.5, labels = A4path4w1)

A1epsolw1 <- round(sum(as.numeric((epsolw1 <= nw) & (epsolw1 >= tff$W)))/(1004/100), digits = 1)

A2epsolw1 <- round(sum(as.numeric((epsolw1 >= nw)&(epsolw1 >= tff$W)))/(1004/100), digits = 1)

A3epsolw1 <- round(sum(as.numeric((epsolw1 >= nw)&(epsolw1 <= tff$W)))/(1004/100), digits=1)

A4epsolw1 <- round(sum(as.numeric((epsolw1 <= nw)&(epsolw1 <= tff$W)))/(1004/100), digits = 1)


plot(nw, epsolw1, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw1)
text(x=0.8, y=0.5, labels = A2epsolw1)
text(x=0.5, y=0.2, labels = A3epsolw1)
text(x=0.2, y=0.5, labels = A4epsolw1)



#####w2  frequency threshold 5, raw consistency threshold 0.964 #####


ttw2 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.964, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw2

#parsimonious solution
psw2 <- eqmcc(ttw2, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw2

#exclude contradictory assumptions & rows: OP*SM*CM + SP*TP*SM*CM (contradict statement of sufficiency)

ttw2new <- ttw2
suf
rows <- findRows("OP*SM*CM + SP*TP*SM*CM", ttw2new, remainders = FALSE)
rows

ttw2new$tt[as.character(rows), "OUT"] <- 0
ttw2new


#conservative solution
csw2 <- eqmcc(ttw2new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw2

#enhanced parsimonious solution (all simplifying assumptions)

epsw2 <- eqmcc(ttw2new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw2

#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  sm*CM + op*SM*cm + sp*op*SM

epsolw2 <- compute("sm*CM + op*SM*cm + sp*op*SM ", data=tff)

obsepsolw2 <- as.numeric(epsolw2 <= nw)
zepsolw2 <- round(((((sum(obsepsolw2 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolw2
pzepsolw2 <- 2*pnorm(-abs(zepsolw2))
pzepsolw2

#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw2 <- eqmcc(ttw2new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw2


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M2: sm*CM + op*SM*cm + sp*op*SM

path1w2 <- compute("sm*CM ", data=tff)
path2w2 <- compute("op*SM*cm", data=tff)
path3w2 <- compute("sp*op*SM ", data=tff)

A1path1w2 <- round(sum(as.numeric((path1w2 <= nw) & (path1w2 >= tff$W)))/(1004/100), digits = 1)

A2path1w2 <- round(sum(as.numeric((path1w2 >= nw)&(path1w2 >= tff$W)))/(1004/100), digits = 1)

A3path1w2 <- round(sum(as.numeric((path1w2 >= nw)&(path1w2 <= tff$W)))/(1004/100), digits=1)

A4path1w2 <- round(sum(as.numeric((path1w2 <= nw)&(path1w2 <= tff$W)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w2)
text(x=0.8, y=0.5, labels = A2path1w2)
text(x=0.5, y=0.2, labels = A3path1w2)
text(x=0.2, y=0.5, labels = A4path1w2)



A1path2w2 <- round(sum(as.numeric((path2w2 <= nw) & (path2w2 >= tff$W)))/(1004/100), digits = 1)

A2path2w2 <- round(sum(as.numeric((path2w2 >= nw)&(path2w2 >= tff$W)))/(1004/100), digits = 1)

A3path2w2 <- round(sum(as.numeric((path2w2 >= nw)&(path2w2 <= tff$W)))/(1004/100), digits=1)

A4path2w2 <- round(sum(as.numeric((path2w2 <= nw)&(path2w2 <= tff$W)))/(1004/100), digits = 1)


plot(nw, path2w2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w2)
text(x=0.8, y=0.5, labels = A2path2w2)
text(x=0.5, y=0.2, labels = A3path2w2)
text(x=0.2, y=0.5, labels = A4path2w2)

A1path3w2 <- round(sum(as.numeric((path3w2 <= nw) & (path3w2 >= tff$W)))/(1004/100), digits = 1)

A2path3w2 <- round(sum(as.numeric((path3w2 >= nw)&(path3w2 >= tff$W)))/(1004/100), digits = 1)

A3path3w2 <- round(sum(as.numeric((path3w2 >= nw)&(path3w2 <= tff$W)))/(1004/100), digits=1)

A4path3w2 <- round(sum(as.numeric((path3w2 <= nw)&(path3w2 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path3w2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 3',
     xlab= 'Path 3',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3w2)
text(x=0.8, y=0.5, labels = A2path3w2)
text(x=0.5, y=0.2, labels = A3path3w2)
text(x=0.2, y=0.5, labels = A4path3w2)


A1epsolw2 <- round(sum(as.numeric((epsolw2 <= nw) & (epsolw2 >= tff$W)))/(1004/100), digits = 1)

A2epsolw2 <- round(sum(as.numeric((epsolw2 >= nw)&(epsolw2 >= tff$W)))/(1004/100), digits = 1)

A3epsolw2 <- round(sum(as.numeric((epsolw2 >= nw)&(epsolw2 <= tff$W)))/(1004/100), digits=1)

A4epsolw2 <- round(sum(as.numeric((epsolw2 <= nw)&(epsolw2 <= tff$W)))/(1004/100), digits = 1)


plot(nw, epsolw2, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw2)
text(x=0.8, y=0.5, labels = A2epsolw2)
text(x=0.5, y=0.2, labels = A3epsolw2)
text(x=0.2, y=0.5, labels = A4epsolw2)

#####w3  frequency threshold 5, raw consistency threshold 0.9613 #####


ttw3 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.9613, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw3

#parsimonious solution
psw3 <- eqmcc(ttw3, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw3

#exclude contradictory assumptions & rows: OP*SM*CM + SP*TP*SM*CM (contradict statement of sufficiency)

ttw3new <- ttw3
suf
rows <- findRows("OP*SM*CM + SP*TP*SM*CM", ttw3new, remainders = FALSE)
rows

ttw3new$tt[as.character(rows), "OUT"] <- 0
ttw3new


#conservative solution
csw3 <- eqmcc(ttw3new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw3

#enhanced parsimonious solution (all simplifying assumptions)

epsw3 <- eqmcc(ttw3new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw3


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  sm*CM + op*SM*cm + sp*TP*op

epsolw3 <- compute("sm*CM + op*SM*cm + sp*TP*op ", data=tff)

obsepsolw3 <- as.numeric(epsolw3 <= nw)
zepsolw3 <- round(((((sum(obsepsolw3 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolw3
pzepsolw3 <- 2*pnorm(-abs(zepsolw3))
pzepsolw3


#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw3 <- eqmcc(ttw3new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw3

# identify simplifying assumptions and easy counterfactuals

SAw3 <- psw3$SA
SAw3

ECw3 <- isw3$i.sol$C1P1$EC
ECw3


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: sm*CM + op*SM*cm + sp*TP*op
path1w3 <- compute("sm*CM", data=tff)
path2w3 <- compute("op*SM*cm ", data=tff)
path3w3 <- compute("sp*TP*op ", data=tff)


A1path1w3 <- round(sum(as.numeric((path1w3 <= nw) & (path1w3 >= tff$W)))/(1004/100), digits = 1)

A2path1w3 <- round(sum(as.numeric((path1w3 >= nw)&(path1w3 >= tff$W)))/(1004/100), digits = 1)

A3path1w3 <- round(sum(as.numeric((path1w3 >= nw)&(path1w3 <= tff$W)))/(1004/100), digits=1)

A4path1w3 <- round(sum(as.numeric((path1w3 <= nw)&(path1w3 <= tff$W)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w3)
text(x=0.8, y=0.5, labels = A2path1w3)
text(x=0.5, y=0.2, labels = A3path1w3)
text(x=0.2, y=0.5, labels = A4path1w3)



A1path2w3 <- round(sum(as.numeric((path2w3 <= nw) & (path2w3 >= tff$W)))/(1004/100), digits = 1)

A2path2w3 <- round(sum(as.numeric((path2w3 >= nw)&(path2w3 >= tff$W)))/(1004/100), digits = 1)

A3path2w3 <- round(sum(as.numeric((path2w3 >= nw)&(path2w3 <= tff$W)))/(1004/100), digits=1)

A4path2w3 <- round(sum(as.numeric((path2w3 <= nw)&(path2w3 <= tff$W)))/(1004/100), digits = 1)


plot(nw, path2w3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w3)
text(x=0.8, y=0.5, labels = A2path2w3)
text(x=0.5, y=0.2, labels = A3path2w3)
text(x=0.2, y=0.5, labels = A4path2w3)

A1path3w3 <- round(sum(as.numeric((path3w3 <= nw) & (path3w3 >= tff$W)))/(1004/100), digits = 1)

A2path3w3 <- round(sum(as.numeric((path3w3 >= nw)&(path3w3 >= tff$W)))/(1004/100), digits = 1)

A3path3w3 <- round(sum(as.numeric((path3w3 >= nw)&(path3w3 <= tff$W)))/(1004/100), digits=1)

A4path3w3 <- round(sum(as.numeric((path3w3 <= nw)&(path3w3 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path3w3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 3',
     xlab= 'Path 3',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3w3)
text(x=0.8, y=0.5, labels = A2path3w3)
text(x=0.5, y=0.2, labels = A3path3w3)
text(x=0.2, y=0.5, labels = A4path3w3)




A1epsolw3 <- round(sum(as.numeric((epsolw3 <= nw) & (epsolw3 >= tff$W)))/(1004/100), digits = 1)

A2epsolw3 <- round(sum(as.numeric((epsolw3 >= nw)&(epsolw3 >= tff$W)))/(1004/100), digits = 1)

A3epsolw3 <- round(sum(as.numeric((epsolw3 >= nw)&(epsolw3 <= tff$W)))/(1004/100), digits=1)

A4epsolw3 <- round(sum(as.numeric((epsolw3 <= nw)&(epsolw3 <= tff$W)))/(1004/100), digits = 1)


plot(nw, epsolw3, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw3)
text(x=0.8, y=0.5, labels = A2epsolw3)
text(x=0.5, y=0.2, labels = A3epsolw3)
text(x=0.2, y=0.5, labels = A4epsolw3)

#Z-Test for intermediate solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  sm*CM + op*SM*cm + sp*TP*op

isolw3 <- compute("sm*CM + op*SM*cm + sp*TP*op ", data=tff)

obsisolw3 <- as.numeric(isolw3 <= nw)
zisolw3 <- round(((((sum(obsisolw3 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zisolw3
pzisolw3 <- 2*pnorm(-abs(zisolw3))
pzisolw3

# for path 1
# Calculate membership in path  sm*CM
path1isolw3 <- compute("sm*CM ", data=tff)

obspath1isolw3 <- as.numeric(path1isolw3 <= nw)
zpath1isolw3 <- round(((((sum(obspath1isolw3 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zpath1isolw3
pzpath1isolw3 <- 2*pnorm(-abs(zpath1isolw3))
pzpath1isolw3

#for path 2
# Calculate membership in path  op*SM*cm
path2isolw3 <- compute("op*SM*cm ", data=tff)

obspath2isolw3 <- as.numeric(path2isolw3 <= nw)
zpath2isolw3 <- round(((((sum(obspath2isolw3 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zpath2isolw3
pzpath2isolw3 <- 2*pnorm(-abs(zpath2isolw3))
pzpath2isolw3

# for path 3
# Calculate membership in path  sp*TP*op
path3isolw3 <- compute("sp*TP*op ", data=tff)

obspath3isolw3 <- as.numeric(path3isolw3 <= nw)
zpath3isolw3 <- round(((((sum(obspath3isolw3 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zpath3isolw3
pzpath3isolw3 <- 2*pnorm(-abs(zpath3isolw3))
pzpath3isolw3


#####w4  frequency threshold 10, raw consistency threshold 0.9611 #####


ttw4 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.9611, n.cut=10, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw4

#parsimonious solution
psw4 <- eqmcc(ttw4, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw4

#exclude contradictory assumptions & rows: OP*SM*CM + SP*TP*SM*CM (contradict statement of sufficiency)

ttw4new <- ttw4
suf
rows <- findRows("OP*SM*CM + SP*TP*SM*CM", ttw4new, remainders = FALSE)
rows

ttw4new$tt[as.character(rows), "OUT"] <- 0
ttw4new


#conservative solution
csw4 <- eqmcc(ttw4new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw4

#enhanced parsimonious solution (all simplifying assumptions)

epsw4 <- eqmcc(ttw4new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw4

#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  sm*CM + op*SM*cm + sp*TP*op

epsolw4 <- compute("sm*CM + op*SM*cm + sp*TP*op ", data=tff)

obsepsolw4 <- as.numeric(epsolw4 <= nw)
zepsolw4 <- round(((((sum(obsepsolw4 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolw4
pzepsolw4 <- 2*pnorm(-abs(zepsolw4))
pzepsolw4


#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw4 <- eqmcc(ttw4new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw4

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: sm*CM + op*SM*cm + sp*TP*op
path1w4 <- compute("sm*CM", data=tff)
path2w4 <- compute("op*SM*cm ", data=tff)
path3w4 <- compute("sp*TP*op ", data=tff)


A1path1w4 <- round(sum(as.numeric((path1w4 <= nw) & (path1w4 >= tff$W)))/(1004/100), digits = 1)

A2path1w4 <- round(sum(as.numeric((path1w4 >= nw)&(path1w4 >= tff$W)))/(1004/100), digits = 1)

A3path1w4 <- round(sum(as.numeric((path1w4 >= nw)&(path1w4 <= tff$W)))/(1004/100), digits=1)

A4path1w4 <- round(sum(as.numeric((path1w4 <= nw)&(path1w4 <= tff$W)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w4)
text(x=0.8, y=0.5, labels = A2path1w4)
text(x=0.5, y=0.2, labels = A3path1w4)
text(x=0.2, y=0.5, labels = A4path1w4)



A1path2w4 <- round(sum(as.numeric((path2w4 <= nw) & (path2w4 >= tff$W)))/(1004/100), digits = 1)

A2path2w4 <- round(sum(as.numeric((path2w4 >= nw)&(path2w4 >= tff$W)))/(1004/100), digits = 1)

A3path2w4 <- round(sum(as.numeric((path2w4 >= nw)&(path2w4 <= tff$W)))/(1004/100), digits=1)

A4path2w4 <- round(sum(as.numeric((path2w4 <= nw)&(path2w4 <= tff$W)))/(1004/100), digits = 1)


plot(nw, path2w4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w4)
text(x=0.8, y=0.5, labels = A2path2w4)
text(x=0.5, y=0.2, labels = A3path2w4)
text(x=0.2, y=0.5, labels = A4path2w4)

A1path3w4 <- round(sum(as.numeric((path3w4 <= nw) & (path3w4 >= tff$W)))/(1004/100), digits = 1)

A2path3w4 <- round(sum(as.numeric((path3w4 >= nw)&(path3w4 >= tff$W)))/(1004/100), digits = 1)

A3path3w4 <- round(sum(as.numeric((path3w4 >= nw)&(path3w4 <= tff$W)))/(1004/100), digits=1)

A4path3w4 <- round(sum(as.numeric((path3w4 <= nw)&(path3w4 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path3w4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 3',
     xlab= 'Path 3',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3w4)
text(x=0.8, y=0.5, labels = A2path3w4)
text(x=0.5, y=0.2, labels = A3path3w4)
text(x=0.2, y=0.5, labels = A4path3w4)

A1path4w4 <- round(sum(as.numeric((path4w4 <= nw) & (path4w4 >= tff$W)))/(1004/100), digits = 1)



A1epsolw4 <- round(sum(as.numeric((epsolw4 <= nw) & (epsolw4 >= tff$W)))/(1004/100), digits = 1)

A2epsolw4 <- round(sum(as.numeric((epsolw4 >= nw)&(epsolw4 >= tff$W)))/(1004/100), digits = 1)

A3epsolw4 <- round(sum(as.numeric((epsolw4 >= nw)&(epsolw4 <= tff$W)))/(1004/100), digits=1)

A4epsolw4 <- round(sum(as.numeric((epsolw4 <= nw)&(epsolw4 <= tff$W)))/(1004/100), digits = 1)


plot(nw, epsolw4, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw4)
text(x=0.8, y=0.5, labels = A2epsolw4)
text(x=0.5, y=0.2, labels = A3epsolw4)
text(x=0.2, y=0.5, labels = A4epsolw4)


#####w5  frequency threshold 10, raw consistency threshold 0.9583 #####


ttw5 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.9583, n.cut=10, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw5

#parsimonious solution
psw5 <- eqmcc(ttw5, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw5

#exclude contradictory assumptions & rows: OOP*SM*CM + SP*TP*SM*CM (contradict statement of sufficiency)

ttw5new <- ttw5
suf
rows <- findRows("OP*SM*CM + SP*TP*SM*CM", ttw5new, remainders = FALSE)
rows

ttw5new$tt[as.character(rows), "OUT"] <- 0
ttw5new


#conservative solution
csw5 <- eqmcc(ttw5new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw5

#enhanced parsimonious solution (all simplifying assumptions)

epsw5 <- eqmcc(ttw5new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw5



#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  sm*CM + op*SM*cm + TP*op*cm 
epsolw5 <- compute("sm*CM + op*SM*cm + TP*op*cm ", data=tff)

obsepsolw5 <- as.numeric(epsolw5 <= nw)
zepsolw5 <- round(((((sum(obsepsolw5 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolw5
pzepsolw5 <- 2*pnorm(-abs(zepsolw5))
pzepsolw5

#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw5 <- eqmcc(ttw5new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw5

#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: sm*CM + op*SM*cm + TP*op*cm
path1w5 <- compute("sm*CM ", data=tff)
path2w5 <- compute("op*SM*cm ", data=tff)
path3w5 <- compute("TP*op*cm ", data=tff)

A1path1w5 <- round(sum(as.numeric((path1w5 <= nw) & (path1w5 >= tff$W)))/(1004/100), digits = 1)

A2path1w5 <- round(sum(as.numeric((path1w5 >= nw)&(path1w5 >= tff$W)))/(1004/100), digits = 1)

A3path1w5 <- round(sum(as.numeric((path1w5 >= nw)&(path1w5 <= tff$W)))/(1004/100), digits=1)

A4path1w5 <- round(sum(as.numeric((path1w5 <= nw)&(path1w5 <= tff$W)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w5)
text(x=0.8, y=0.5, labels = A2path1w5)
text(x=0.5, y=0.2, labels = A3path1w5)
text(x=0.2, y=0.5, labels = A4path1w5)

A1path2w5 <- round(sum(as.numeric((path2w5 <= nw) & (path2w5 >= tff$W)))/(1004/100), digits = 1)

A2path2w5 <- round(sum(as.numeric((path2w5 >= nw)&(path2w5 >= tff$W)))/(1004/100), digits = 1)

A3path2w5 <- round(sum(as.numeric((path2w5 >= nw)&(path2w5 <= tff$W)))/(1004/100), digits=1)

A4path2w5 <- round(sum(as.numeric((path2w5 <= nw)&(path2w5 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path2w5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w5)
text(x=0.8, y=0.5, labels = A2path2w5)
text(x=0.5, y=0.2, labels = A3path2w5)
text(x=0.2, y=0.5, labels = A4path2w5)

A1path3w5 <- round(sum(as.numeric((path3w5 <= nw) & (path3w5 >= tff$W)))/(1004/100), digits = 1)

A2path3w5 <- round(sum(as.numeric((path3w5 >= nw)&(path3w5 >= tff$W)))/(1004/100), digits = 1)

A3path3w5 <- round(sum(as.numeric((path3w5 >= nw)&(path3w5 <= tff$W)))/(1004/100), digits=1)

A4path3w5 <- round(sum(as.numeric((path3w5 <= nw)&(path3w5 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path3w5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 3',
     xlab= 'Path 3',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3w5)
text(x=0.8, y=0.5, labels = A2path3w5)
text(x=0.5, y=0.2, labels = A3path3w5)
text(x=0.2, y=0.5, labels = A4path3w5)


A1epsolw5 <- round(sum(as.numeric((epsolw5 <= nw) & (epsolw5 >= tff$W)))/(1004/100), digits = 1)

A2epsolw5 <- round(sum(as.numeric((epsolw5 >= nw)&(epsolw5 >= tff$W)))/(1004/100), digits = 1)

A3epsolw5 <- round(sum(as.numeric((epsolw5 >= nw)&(epsolw5 <= tff$W)))/(1004/100), digits=1)

A4epsolw5 <- round(sum(as.numeric((epsolw5 <= nw)&(epsolw5 <= tff$W)))/(1004/100), digits = 1)


plot(nw, epsolw5, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw5)
text(x=0.8, y=0.5, labels = A2epsolw5)
text(x=0.5, y=0.2, labels = A3epsolw5)
text(x=0.2, y=0.5, labels = A4epsolw5)


#####w6  frequency threshold 10, raw consistency threshold 0.956#####


ttw6 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.956, n.cut=10, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw6

#parsimonious solution
psw6 <- eqmcc(ttw6, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw6

#exclude contradictory assumptions & rows: OP*SM*CM + SP*TP*SM*CM  (contradict statement of sufficiency)

ttw6new <- ttw6
suf
rows <- findRows("OP*SM*CM + SP*TP*SM*CM ", ttw6new, remainders = FALSE)
rows

ttw6new$tt[as.character(rows), "OUT"] <- 0
ttw6new


#conservative solution
csw6 <- eqmcc(ttw6new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw6

#enhanced parsimonious solution (all simplifying assumptions)

epsw6 <- eqmcc(ttw6new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw6



#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  sm*CM + op*SM*cm + TP*op*cm + sp*TP*SM*cm
epsolw6 <- compute("sm*CM + op*SM*cm + TP*op*cm + sp*TP*SM*cm ", data=tff)

obsepsolw6 <- as.numeric(epsolw6 <= nw)
zepsolw6 <- round(((((sum(obsepsolw6 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolw6
pzepsolw6 <- 2*pnorm(-abs(zepsolw6))
pzepsolw6


#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw6 <- eqmcc(ttw6new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw6



#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: sm*CM + op*SM*cm + TP*op*cm + sp*TP*SM*cm
path1w6 <- compute("sm*CM ", data=tff)
path2w6 <- compute("op*SM*cm ", data=tff)
path3w6 <- compute("TP*op*cm ", data=tff)
path4w6 <- compute("sp*TP*SM*cm ", data=tff)

A1path1w6 <- round(sum(as.numeric((path1w6 <= nw) & (path1w6 >= tff$W)))/(1004/100), digits = 1)

A2path1w6 <- round(sum(as.numeric((path1w6 >= nw)&(path1w6 >= tff$W)))/(1004/100), digits = 1)

A3path1w6 <- round(sum(as.numeric((path1w6 >= nw)&(path1w6 <= tff$W)))/(1004/100), digits=1)

A4path1w6 <- round(sum(as.numeric((path1w6 <= nw)&(path1w6 <= tff$W)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w6)
text(x=0.8, y=0.5, labels = A2path1w6)
text(x=0.5, y=0.2, labels = A3path1w6)
text(x=0.2, y=0.5, labels = A4path1w6)

A1path2w6 <- round(sum(as.numeric((path2w6 <= nw) & (path2w6 >= tff$W)))/(1004/100), digits = 1)

A2path2w6 <- round(sum(as.numeric((path2w6 >= nw)&(path2w6 >= tff$W)))/(1004/100), digits = 1)

A3path2w6 <- round(sum(as.numeric((path2w6 >= nw)&(path2w6 <= tff$W)))/(1004/100), digits=1)

A4path2w6 <- round(sum(as.numeric((path2w6 <= nw)&(path2w6 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path2w6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w6)
text(x=0.8, y=0.5, labels = A2path2w6)
text(x=0.5, y=0.2, labels = A3path2w6)
text(x=0.2, y=0.5, labels = A4path2w6)

A1path3w6 <- round(sum(as.numeric((path3w6 <= nw) & (path3w6 >= tff$W)))/(1004/100), digits = 1)

A2path3w6 <- round(sum(as.numeric((path3w6 >= nw)&(path3w6 >= tff$W)))/(1004/100), digits = 1)

A3path3w6 <- round(sum(as.numeric((path3w6 >= nw)&(path3w6 <= tff$W)))/(1004/100), digits=1)

A4path3w6 <- round(sum(as.numeric((path3w6 <= nw)&(path3w6 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path3w6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 3',
     xlab= 'Path 3',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3w6)
text(x=0.8, y=0.5, labels = A2path3w6)
text(x=0.5, y=0.2, labels = A3path3w6)
text(x=0.2, y=0.5, labels = A4path3w6)


A1path4w6 <- round(sum(as.numeric((path4w6 <= nw) & (path4w6 >= tff$W)))/(1004/100), digits = 1)

A2path4w6 <- round(sum(as.numeric((path4w6 >= nw)&(path4w6 >= tff$W)))/(1004/100), digits = 1)

A3path4w6 <- round(sum(as.numeric((path4w6 >= nw)&(path4w6 <= tff$W)))/(1004/100), digits=1)

A4path4w6 <- round(sum(as.numeric((path4w6 <= nw)&(path4w6 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path4w6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 4',
     xlab= 'Path 4',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path4w6)
text(x=0.8, y=0.5, labels = A2path4w6)
text(x=0.5, y=0.2, labels = A3path4w6)
text(x=0.2, y=0.5, labels = A4path4w6)


A1epsolw6 <- round(sum(as.numeric((epsolw6 <= nw) & (epsolw6 >= tff$W)))/(1004/100), digits = 1)

A2epsolw6 <- round(sum(as.numeric((epsolw6 >= nw)&(epsolw6 >= tff$W)))/(1004/100), digits = 1)

A3epsolw6 <- round(sum(as.numeric((epsolw6 >= nw)&(epsolw6 <= tff$W)))/(1004/100), digits=1)

A4epsolw6 <- round(sum(as.numeric((epsolw6 <= nw)&(epsolw6 <= tff$W)))/(1004/100), digits = 1)


plot(nw, epsolw6, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw6)
text(x=0.8, y=0.5, labels = A2epsolw6)
text(x=0.5, y=0.2, labels = A3epsolw6)
text(x=0.2, y=0.5, labels = A4epsolw6)

#####w7  frequency threshold 50, raw consistency threshold 0.918#####


ttw7 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.918, n.cut=50, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw7

#parsimonious solution
psw7 <- eqmcc(ttw7, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw7

#exclude contradictory assumptions & rows: OP*SM*CM + SP*TP*SM*CM (contradict statement of sufficiency)

ttw7new <- ttw7
suf
rows <- findRows("OP*SM*CM + SP*TP*SM*CM", ttw7new, remainders = FALSE)
rows

ttw7new$tt[as.character(rows), "OUT"] <- 0
ttw7new


#conservative solution
csw7 <- eqmcc(ttw7new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw7

#enhanced parsimonious solution (all simplifying assumptions)

epsw7 <- eqmcc(ttw7new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw7


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  op*sm + sp*TP*sm + SP*tp*sm
epsolw7 <- compute("op*sm + sp*TP*sm + SP*tp*sm", data=tff)

obsepsolw7 <- as.numeric(epsolw7 <= nw)
zepsolw7 <- round(((((sum(obsepsolw7 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolw7
pzepsolw7 <- 2*pnorm(-abs(zepsolw7))
pzepsolw7

#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw7 <- eqmcc(ttw7new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw7


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M8: op*sm + sp*TP*sm + SP*tp*sm
path1w7 <- compute("op*sm ", data=tff)
path2w7 <- compute("sp*TP*sm ", data=tff)
path3w7 <- compute("SP*tp*sm ", data=tff)

A1path1w7 <- round(sum(as.numeric((path1w7 <= nw) & (path1w7 >= tff$W)))/(1004/100), digits = 1)

A2path1w7 <- round(sum(as.numeric((path1w7 >= nw)&(path1w7 >= tff$W)))/(1004/100), digits = 1)

A3path1w7 <- round(sum(as.numeric((path1w7 >= nw)&(path1w7 <= tff$W)))/(1004/100), digits=1)

A4path1w7 <- round(sum(as.numeric((path1w7 <= nw)&(path1w7 <= tff$W)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w7, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w7)
text(x=0.8, y=0.5, labels = A2path1w7)
text(x=0.5, y=0.2, labels = A3path1w7)
text(x=0.2, y=0.5, labels = A4path1w7)

A1path2w7 <- round(sum(as.numeric((path2w7 <= nw) & (path2w7 >= tff$W)))/(1004/100), digits = 1)

A2path2w7 <- round(sum(as.numeric((path2w7 >= nw)&(path2w7 >= tff$W)))/(1004/100), digits = 1)

A3path2w7 <- round(sum(as.numeric((path2w7 >= nw)&(path2w7 <= tff$W)))/(1004/100), digits=1)

A4path2w7 <- round(sum(as.numeric((path2w7 <= nw)&(path2w7 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path2w7, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w7)
text(x=0.8, y=0.5, labels = A2path2w7)
text(x=0.5, y=0.2, labels = A3path2w7)
text(x=0.2, y=0.5, labels = A4path2w7)

A1path3w7 <- round(sum(as.numeric((path3w7 <= nw) & (path3w7 >= tff$W)))/(1004/100), digits = 1)

A2path3w7 <- round(sum(as.numeric((path3w7 >= nw)&(path3w7 >= tff$W)))/(1004/100), digits = 1)

A3path3w7 <- round(sum(as.numeric((path3w7 >= nw)&(path3w7 <= tff$W)))/(1004/100), digits=1)

A4path3w7 <- round(sum(as.numeric((path3w7 <= nw)&(path3w7 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path3w7, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 3',
     xlab= 'Path 3',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path3w7)
text(x=0.8, y=0.5, labels = A2path3w7)
text(x=0.5, y=0.2, labels = A3path3w7)
text(x=0.2, y=0.5, labels = A4path3w7)


A1epsolw7 <- round(sum(as.numeric((epsolw7 <= nw) & (epsolw7 >= tff$W)))/(1004/100), digits = 1)

A2epsolw7 <- round(sum(as.numeric((epsolw7 >= nw)&(epsolw7 >= tff$W)))/(1004/100), digits = 1)

A3epsolw7 <- round(sum(as.numeric((epsolw7 >= nw)&(epsolw7 <= tff$W)))/(1004/100), digits=1)

A4epsolw7 <- round(sum(as.numeric((epsolw7 <= nw)&(epsolw7 <= tff$W)))/(1004/100), digits = 1)


plot(nw, epsolw7, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw7)
text(x=0.8, y=0.5, labels = A2epsolw7)
text(x=0.5, y=0.2, labels = A3epsolw7)
text(x=0.2, y=0.5, labels = A4epsolw7)

#####w8  frequency threshold 50, raw consistency threshold 0.901 #####


ttw8 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.901, n.cut=50, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw8

#parsimonious solution
psw8 <- eqmcc(ttw8, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw8

#exclude contradictory assumptions & rows: OP*SM*CM + SP*TP*SM*CM (contradict statement of sufficiency)

ttw8new <- ttw8
suf
rows <- findRows("OP*SM*CM + SP*TP*SM*CM", ttw8new, remainders = FALSE)
rows

ttw8new$tt[as.character(rows), "OUT"] <- 0
ttw8new


#conservative solution
csw8 <- eqmcc(ttw8new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw8

#enhanced parsimonious solution (all simplifying assumptions)

epsw8 <- eqmcc(ttw8new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw8


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  sp*sm + tp*sm
epsolw8 <- compute("sp*sm + tp*sm ", data=tff)

obsepsolw8 <- as.numeric(epsolw8 <= nw)
zepsolw8 <- round(((((sum(obsepsolw8 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolw8
pzepsolw8 <- 2*pnorm(-abs(zepsolw8))
pzepsolw8

#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw8 <- eqmcc(ttw8new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw8


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M4: sp*sm + tp*sm
path1w8 <- compute("sp*sm ", data=tff)
path2w8 <- compute("tp*sm ", data=tff)

A1path1w8 <- round(sum(as.numeric((path1w8 <= nw) & (path1w8 >= tff$W)))/(1004/100), digits = 1)

A2path1w8 <- round(sum(as.numeric((path1w8 >= nw)&(path1w8 >= tff$W)))/(1004/100), digits = 1)

A3path1w8 <- round(sum(as.numeric((path1w8 >= nw)&(path1w8 <= tff$W)))/(1004/100), digits=1)

A4path1w8 <- round(sum(as.numeric((path1w8 <= nw)&(path1w8 <= tff$W)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w8, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w8)
text(x=0.8, y=0.5, labels = A2path1w8)
text(x=0.5, y=0.2, labels = A3path1w8)
text(x=0.2, y=0.5, labels = A4path1w8)

A1path2w8 <- round(sum(as.numeric((path2w8 <= nw) & (path2w8 >= tff$W)))/(1004/100), digits = 1)

A2path2w8 <- round(sum(as.numeric((path2w8 >= nw)&(path2w8 >= tff$W)))/(1004/100), digits = 1)

A3path2w8 <- round(sum(as.numeric((path2w8 >= nw)&(path2w8 <= tff$W)))/(1004/100), digits=1)

A4path2w8 <- round(sum(as.numeric((path2w8 <= nw)&(path2w8 <= tff$W)))/(1004/100), digits = 1)

plot(nw, path2w8, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 2',
     xlab= 'Path 2',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path2w8)
text(x=0.8, y=0.5, labels = A2path2w8)
text(x=0.5, y=0.2, labels = A3path2w8)
text(x=0.2, y=0.5, labels = A4path2w8)

A1epsolw8 <- round(sum(as.numeric((epsolw8 <= nw) & (epsolw8 >= tff$W)))/(1004/100), digits = 1)

A2epsolw8 <- round(sum(as.numeric((epsolw8 >= nw)&(epsolw8 >= tff$W)))/(1004/100), digits = 1)

A3epsolw8 <- round(sum(as.numeric((epsolw8 >= nw)&(epsolw8 <= tff$W)))/(1004/100), digits=1)

A4epsolw8 <- round(sum(as.numeric((epsolw8 <= nw)&(epsolw8 <= tff$W)))/(1004/100), digits = 1)


plot(nw, epsolw8, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw8)
text(x=0.8, y=0.5, labels = A2epsolw8)
text(x=0.5, y=0.2, labels = A3epsolw8)
text(x=0.2, y=0.5, labels = A4epsolw8)


#####w9  frequency threshold 50, raw consistency threshold 0.899#####


ttw9 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.899, n.cut=50, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw9

#parsimonious solution
psw9 <- eqmcc(ttw9, include="?", details=TRUE, show.cases=FALSE,  row.dom=TRUE, all.sol=FALSE)

psw9

#exclude contradictory assumptions & rows: OP*SM*CM + SP*TP*SM*CM (contradict statement of sufficiency)

ttw9new <- ttw9
suf
rows <- findRows("OP*SM*CM + SP*TP*SM*CM", ttw9new, remainders = FALSE)
rows

ttw9new$tt[as.character(rows), "OUT"] <- 0
ttw9new


#conservative solution
csw9 <- eqmcc(ttw9new, details=TRUE, PRI=TRUE, show.cases=FALSE, rowdom=TRUE, use.tilde=FALSE)
csw9

#enhanced parsimonious solution (all simplifying assumptions)

epsw9 <- eqmcc(ttw9new, include="?", details=TRUE, PRI=TRUE, show.cases=FALSE, 
               rowdom=TRUE, use.tilde=FALSE)
epsw9


#Z-Test for enhanced parsimonious solution, 0.8 (almost always sufficient)
# Calculate membership in solution term  sm
epsolw9 <- compute("sm ", data=tff)

obsepsolw9 <- as.numeric(epsolw9 <= nw)
zepsolw9 <- round(((((sum(obsepsolw9 == 1))/1004) - 0.8)-(1/(2*1004))) / (sqrt((0.8*(1-0.8))/1004)), digits=1)
zepsolw9
pzepsolw9 <- 2*pnorm(-abs(zepsolw9))
pzepsolw9

#intermediate solution with Directional expectations: sp->w, tp->w, op->w.
isw9 <- eqmcc(ttw9new, details = TRUE, include = "?",  use.tilde = FALSE, PRI = TRUE,
              row.dom = TRUE, show.cases = FALSE, 
              dir.exp = "0, 0, 0, -, -")

isw9


#Test for skewness of sufficient conditions (enhanced parsimonious solution)
# calculate membership in paths of M1: sm
path1w9 <- compute("sm ", data=tff)


A1path1w9 <- round(sum(as.numeric((path1w9 <= nw) & (path1w9 >= tff$W)))/(1004/100), digits = 1)

A2path1w9 <- round(sum(as.numeric((path1w9 >= nw)&(path1w9 >= tff$W)))/(1004/100), digits = 1)

A3path1w9 <- round(sum(as.numeric((path1w9 >= nw)&(path1w9 <= tff$W)))/(1004/100), digits=1)

A4path1w9 <- round(sum(as.numeric((path1w9 <= nw)&(path1w9 <= tff$W)))/(1004/100), digits = 1)


par(mfrow=c(2, 3))
plot(nw, path1w9, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness path 1',
     xlab= 'Path 1',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1path1w9)
text(x=0.8, y=0.5, labels = A2path1w9)
text(x=0.5, y=0.2, labels = A3path1w9)
text(x=0.2, y=0.5, labels = A4path1w9)


A1epsolw9 <- round(sum(as.numeric((epsolw9 <= nw) & (epsolw9 >= tff$W)))/(1004/100), digits = 1)

A2epsolw9 <- round(sum(as.numeric((epsolw9 >= nw)&(epsolw9 >= tff$W)))/(1004/100), digits = 1)

A3epsolw9 <- round(sum(as.numeric((epsolw9 >= nw)&(epsolw9 <= tff$W)))/(1004/100), digits=1)

A4epsolw9 <- round(sum(as.numeric((epsolw9 <= nw)&(epsolw9 <= tff$W)))/(1004/100), digits = 1)


plot(nw, epsolw9, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main= 'Skewness eps',
     xlab= 'solution term',
     ylab='w')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1epsolw9)
text(x=0.8, y=0.5, labels = A2epsolw9)
text(x=0.5, y=0.2, labels = A3epsolw9)
text(x=0.2, y=0.5, labels = A4epsolw9)





###############################THEORY EVALUATION COVERAGE ############

tff <- read.csv("tff_d1_c1.csv", row.names=1, header = TRUE, sep = ";",  dec = ",")
describe(tff)
nw <- 1-tff$W



## Theory evaluation for hypothesis 2####

## calculate intersections

Tw <- "sp*tp*op"
Tw
tw <- deMorgan("sp*tp*op")
tw
Sw <- "sm*CM + op*SM*cm + sp*TP*op"
Sw
sw <- deMorgan("sm*CM + op*SM*cm + sp*TP*op")
sw

TSw <- intersection("sp*tp*op", "sm*CM + op*SM*cm + sp*TP*op", snames = "SP, TP, OP, SM, CM")
TSw

tw
tSw <- intersection("OP + SP + TP", "sm*CM + op*SM*cm + sp*TP*op", snames = "SP, TP, OP, SM, CM")
tSw

sw
tsw <- intersection("OP + SP + TP", "cm*OP + cm*sm*SP + cm*sm*tp + CM*SM*SP + CM*SM*tp + OP*SM", snames = "SP, TP, OP, SM, CM")
tsw

Tw
sw
Tsw <- intersection("sp*tp*op", "cm*OP + cm*sm*SP + cm*sm*tp + CM*SM*SP + CM*SM*tp + OP*SM", snames = "SP, TP, OP, SM, CM")
Tsw

#calculate membership in intersections

TSw
TSw1 <- compute("sp*tp*op*sm*CM", data=tff)
TSw2 <- compute("sp*tp*op*SM*cm", data=tff)
TSw <- compute("sp*tp*op*sm*CM + sp*tp*op*SM*cm", data=tff)

Tsw
Tsw1 <- compute("sp*tp*op*sm*cm", data=tff)
Tsw2 <- compute("sp*tp*op*SM*CM", data=tff)
Tsw <- compute("sp*tp*op*sm*cm + sp*tp*op*SM*CM", data=tff)


tSw
tSw1 <- compute("OP*sm*CM", data=tff)
tSw2 <- compute("SP*sm*CM", data=tff)
tSw3 <- compute("SP*op*SM*cm", data=tff)
tSw4 <- compute("TP*sm*CM", data=tff)
tSw5 <- compute("TP*op*SM*cm", data=tff)
tSw6 <- compute("sp*TP*op", data=tff)
tSw <- compute("OP*sm*CM + SP*sm*CM + SP*op*SM*cm + TP*sm*CM + TP*op*SM*cm + sp*TP*op", data=tff)



tsw
tsw1 <- compute("OP*cm", data=tff)
tsw2 <- compute("SP*OP*sm*cm", data=tff)
tsw3 <- compute("tp*OP*sm*cm", data=tff)
tsw4 <- compute("SP*OP*SM*CM", data=tff)
tsw5 <- compute("tp*OP*SM*CM", data=tff)
tsw6 <- compute("OP*SM", data=tff)
tsw7 <- compute("SP*OP*cm", data=tff)
tsw8 <- compute("SP*sm*cm", data=tff)
tsw9 <- compute("SP*tp*sm*cm", data=tff)
tsw10 <- compute("SP*SM*CM", data=tff)
tsw11 <- compute("SP*tp*SM*CM", data=tff)
tsw12 <- compute("SP*OP*SM", data=tff)
tsw13 <- compute("TP*OP*cm", data=tff)
tsw14 <- compute("SP*TP*sm*cm", data=tff)
tsw15 <- compute("SP*TP*SM*CM", data=tff)
tsw16 <- compute("TP*OP*SM", data=tff)
tsw <- compute("OP*cm + SP*OP*sm*cm + tp*OP*sm*cm + SP*OP*SM*CM + tp*OP*SM*CM + OP*SM + SP*OP*cm + SP*sm*cm + SP*tp*sm*cm + SP*SM*CM + SP*tp*SM*CM + SP*OP*SM + TP*OP*cm + SP*TP*sm*cm + SP*TP*SM*CM + TP*OP*SM", data=tff)



#check for logical remainders
TSw
sum(as.numeric(TSw1 > 0.5))
sum(as.numeric(TSw2 > 0.5))
sum(as.numeric(TSw > 0.5))

Tsw
sum(as.numeric(Tsw1 > 0.5))
sum(as.numeric(Tsw2 > 0.5))
sum(as.numeric(Tsw > 0.5))

tSw
sum(as.numeric(tSw1 > 0.5))
sum(as.numeric(tSw2 > 0.5))
sum(as.numeric(tSw3 > 0.5))
sum(as.numeric(tSw4 > 0.5))
sum(as.numeric(tSw5 > 0.5))
sum(as.numeric(tSw6 > 0.5))
sum(as.numeric(tSw > 0.5))


sum(as.numeric(tsw1 > 0.5))
sum(as.numeric(tsw2 > 0.5))
sum(as.numeric(tsw3 > 0.5))
sum(as.numeric(tsw4 > 0.5))
sum(as.numeric(tsw5 > 0.5))
sum(as.numeric(tsw6 > 0.5))
sum(as.numeric(tsw7 > 0.5))
sum(as.numeric(tsw8 > 0.5))
sum(as.numeric(tsw9 > 0.5))
sum(as.numeric(tsw10 > 0.5))
sum(as.numeric(tsw11 > 0.5))
sum(as.numeric(tsw12 > 0.5))
sum(as.numeric(tsw13 > 0.5))
sum(as.numeric(tsw14 > 0.5))
sum(as.numeric(tsw15 > 0.5))
sum(as.numeric(tsw16 > 0.5))
sum(as.numeric(tsw > 0.5))


#calculate % with W and w

TSwW <- round(sum(as.numeric((TSw > 0.5) & (tff$W > 0.5)))/(1004/100), digits = 1)
TSwW

TSww <- round(sum(as.numeric((TSw > 0.5 )&(tff$W < 0.5)))/(1004/100), digits = 1)
TSww

TswW <- round(sum(as.numeric((Tsw > 0.5) & (tff$W > 0.5)))/(1004/100), digits = 1)
TswW

Tsww <- round(sum(as.numeric((Tsw > 0.5 )&(tff$W < 0.5)))/(1004/100), digits = 1)
Tsww

tswW <- round(sum(as.numeric((tsw > 0.5) & (tff$W > 0.5)))/(1004/100), digits = 1)
tswW

tsww <- round(sum(as.numeric((tsw > 0.5 )&(tff$W < 0.5)))/(1004/100), digits = 1)
tsww

tSwW <- round(sum(as.numeric((tSw > 0.5) & (tff$W > 0.5)))/(1004/100), digits = 1)
tSwW

tSww <- round(sum(as.numeric((tSw > 0.5 )&(tff$W < 0.5)))/(1004/100), digits = 1)
tSww



#check for mistakes
fulltew <- sum(TSwW, TSww, tSwW, tSww, TswW, Tsww, tswW, tsww)
fulltew



####Theory evaluation for hypothesis 3#####
TW <- " OP*SM + OP*CM + SP*SM + SP*CM + TP*SM + TP*CM "
TW
tW <- deMorgan("OP*SM + OP*CM + SP*SM + SP*CM + TP*SM + TP*CM ")
tW
SW <- "OP*SM*CM + SP*TP*SM*CM"
SW
sW <- deMorgan("OP*SM*CM + SP*TP*SM*CM")
sW

TW
SW
TSW <- intersection("OP*SM + OP*CM + SP*SM + SP*CM + TP*SM + TP*CM", " OP*SM*CM + SP*TP*SM*CM", snames = "SP, TP, OP, SM, CM")
TSW


tW
SW
tSW <- intersection("op*sp*tp + cm*sm", "OP*SM*CM + SP*TP*SM*CM", snames = "SP, TP, OP, SM, CM")
tSW


sW
tW
tsW <- intersection("op*sp*tp + cm*sm", "cm + op*sp + op*tp + sm", snames = "SP, TP, OP, SM, CM")
tsW


TW
sW
TsW <- intersection("OP*SM + OP*CM + SP*SM + SP*CM + TP*SM + TP*CM", "cm + op*sp + op*tp + sm", snames = "SP, TP, OP, SM, CM")
TsW


#calculate membership in intersections


TSW
TSW1 <- compute("OP*SM*CM", data=tff)
TSW2 <- compute("SP*TP*OP*SM*CM", data=tff)
TSW3 <- compute("SP*OP*SM*CM", data=tff)
TSW4 <- compute("SP*TP*SM*CM", data=tff)
TSW5 <- compute("TP*OP*SM*CM", data=tff)

TSW <- compute("OP*SM*CM + SP*TP*OP*SM*CM + SP*OP*SM*CM + SP*TP*SM*CM + TP*OP*SM*CM", data=tff)

TsW
TsW1 <- compute("OP*SM*cm", data=tff)
TsW2 <- compute("OP*sm*CM", data=tff)
TsW3 <- compute("SP*SM*cm", data=tff)
TsW4 <- compute("SP*tp*op*SM", data=tff)
TsW5 <- compute("SP*tp*op*CM", data=tff)
TsW6 <- compute("SP*sm*CM", data=tff)
TsW7 <- compute("TP*SM*cm", data=tff)
TsW8 <- compute("sp*TP*op*SM", data=tff)
TsW9 <- compute("sp*TP*op*CM", data=tff)
TsW10 <- compute("TP*sm*CM", data=tff)

TsW <- compute("OP*SM*cm + OP*sm*CM + SP*SM*cm + SP*tp*op*SM + SP*tp*op*CM + SP*sm*CM + TP*SM*cm + sp*TP*op*SM + sp*TP*op*CM + TP*sm*CM", data=tff)

tSW

tsW
tsW1 <- compute("sp*tp*op*cm", data=tff)
tsW2 <- compute("sp*tp*op", data=tff)
tsW3 <- compute("sp*tp*op*sm", data=tff)
tsW4 <- compute("sm*cm", data=tff)
tsW5 <- compute("sp*op*sm*cm", data=tff)
tsW6 <- compute("tp*op*sm*cm", data=tff)

tsW <- compute("sp*tp*op*cm + sp*tp*op + sp*tp*op*sm + sm*cm + sp*op*sm*cm + tp*op*sm*cm", data=tff)

#check for logical remainders
TSW
sum(as.numeric(TSW1 > 0.5))
sum(as.numeric(TSW2 > 0.5))
sum(as.numeric(TSW3 > 0.5))
sum(as.numeric(TSW4 > 0.5))
sum(as.numeric(TSW5 > 0.5))
sum(as.numeric(TSW > 0.5))

TsW
sum(as.numeric(TsW1 > 0.5))
sum(as.numeric(TsW2 > 0.5))
sum(as.numeric(TsW3 > 0.5))
sum(as.numeric(TsW4 > 0.5))
sum(as.numeric(TsW5 > 0.5))
sum(as.numeric(TsW6 > 0.5))
sum(as.numeric(TsW7 > 0.5))
sum(as.numeric(TsW8 > 0.5))
sum(as.numeric(TsW9 > 0.5))
sum(as.numeric(TsW10 > 0.5))
sum(as.numeric(TsW > 0.5))

tsW
sum(as.numeric(tsW1 > 0.5))
sum(as.numeric(tsW2 > 0.5))
sum(as.numeric(tsW3 > 0.5))
sum(as.numeric(tsW4 > 0.5))
sum(as.numeric(tsW5 > 0.5))
sum(as.numeric(tsW6 > 0.5))

sum(as.numeric(tsW > 0.5))

#calculate % with W and w

TSWW <- round(sum(as.numeric((TSW > 0.5) & (tff$W > 0.5)))/(1004/100), digits = 1)
TSWW

TSWw <- round(sum(as.numeric((TSW > 0.5 )&(tff$W < 0.5)))/(1004/100), digits = 1)
TSWw

TsWW <- round(sum(as.numeric((TsW > 0.5) & (tff$W > 0.5)))/(1004/100), digits = 1)
TsWW

TsWw <- round(sum(as.numeric((TsW > 0.5 )&(tff$W < 0.5)))/(1004/100), digits = 1)
TsWw

tsWW <- round(sum(as.numeric((tsW > 0.5) & (tff$W > 0.5)))/(1004/100), digits = 1)
tsWW

tsWw <- round(sum(as.numeric((tsW > 0.5 )&(tff$W < 0.5)))/(1004/100), digits = 1)
tsWw

tSWW <- round(sum(as.numeric((tSW > 0.5) & (tff$W > 0.5)))/(1004/100), digits = 1)
tSWW

tSWw <- round(sum(as.numeric((tSW > 0.5 )&(tff$W < 0.5)))/(1004/100), digits = 1)
tSWw



#check for mistakes
fulltew <- sum(TSWW, TSWw, tSWW, tSWw, TsWW, TsWw, tsWW, tsWw)
fulltew



